diff --git a/Exchange.cxx b/Exchange.cxx
index 26cbcb05..644d9ea3 100644
--- a/Exchange.cxx
+++ b/Exchange.cxx
@@ -406,6 +406,7 @@ void cxxExchange::totalize()
   for (std::list<cxxExchComp>::const_iterator it = exchComps.begin(); it != exchComps.end(); ++it) 
   {
     this->totals.add_extensive(it->get_totals(), 1.0);
+    this->totals.add("Charge", it->get_charge_balance());
   }
   return;
 }
diff --git a/ISolution.cxx b/ISolution.cxx
index 93f9e5f2..49e8fa58 100644
--- a/ISolution.cxx
+++ b/ISolution.cxx
@@ -16,7 +16,7 @@
 #include <cassert>     // assert
 #include <algorithm>   // std::sort 
 #include <sstream>
-
+extern void ORCH_write_chemistry_species(std::ostream &chemistry_dat);
 //////////////////////////////////////////////////////////////////////
 // Construction/Destruction
 //////////////////////////////////////////////////////////////////////
@@ -426,7 +426,7 @@ void cxxISolution::ORCH_write_chemistry(std::ostream &chemistry_dat)
   this->ORCH_write_chemistry_primary(chemistry_dat);
   this->ORCH_write_chemistry_total_O_H(chemistry_dat);
   this->ORCH_write_chemistry_alkalinity(chemistry_dat);
-  this->ORCH_write_chemistry_species(chemistry_dat);
+  ORCH_write_chemistry_species(chemistry_dat);
   this->ORCH_write_chemistry_minerals(chemistry_dat);
 }
 
@@ -587,7 +587,6 @@ void cxxISolution::ORCH_write_chemistry_total_O_H(std::ostream &chemistry_dat)
   chemistry_dat << "\")" << std::endl;
   chemistry_dat << std::endl;
 }
-
 void cxxISolution::ORCH_write_chemistry_alkalinity(std::ostream &chemistry_dat)
 {
   chemistry_dat << std::endl << "//********* Alkalinity definitions" << std::endl;
@@ -656,41 +655,7 @@ void cxxISolution::ORCH_write_chemistry_alkalinity(std::ostream &chemistry_dat)
     }
   }
 }
-void cxxISolution::ORCH_write_chemistry_species(std::ostream &chemistry_dat)
-{
-  chemistry_dat << std::endl << "//********* Aqueous species" << std::endl;
-  //
-  //  Write species definitions 
-  //
-  for (int i = 0; i < count_s_x; i++) {
-    if (s_x[i]->primary != NULL) continue;
-    if (s_x[i]->secondary != NULL && s_x[i]->secondary->in == TRUE) continue;
 
-    // write species
-    std::string name(s_x[i]->name);
-    std::replace(name.begin(), name.end(), '(', '[');
-    std::replace(name.begin(), name.end(), ')', ']');
-    chemistry_dat << "@species(" << name.c_str() << ", " << s_x[i]->z << ")" << std::endl;
-
-    // write reaction
-    chemistry_dat << "@reaction(" << name.c_str() << ", " << pow(10.0,  s_x[i]->lk);
-    struct rxn_token *next_token = s_x[i]->rxn_x->token;
-    next_token++;
-    while (next_token->s != NULL || next_token->name != NULL)
-    {
-      chemistry_dat << ", " << next_token->coef;
-      if (next_token->s != NULL)
-      {
-	chemistry_dat << ", " << next_token->s->name;
-      } else {
-	chemistry_dat << ", " << next_token->name;
-      }
-      next_token++;
-    }
-    chemistry_dat << ")" << std::endl;
-  }
-  chemistry_dat << std::endl;
-}
 void cxxISolution::ORCH_write_chemistry_minerals(std::ostream &chemistry_dat)
 {
   chemistry_dat << std::endl << "//********* Adjustments to mineral equilibrium" << std::endl;
@@ -881,12 +846,13 @@ void cxxISolution::ORCH_write_output_vars(std::ostream &outstream)
   //cb
   outstream << "\tchargebalance";
   //mass_water;
-  outstream << "\tH2O.diss";
+  outstream << "\tH2O.con";
   //total_alkalinity;
   outstream << "\tAlkalinity";
   //orchestra master variables
   outstream << "\tH+.diss";
   outstream << "\te-.diss";
+  outstream << "\tH2O.diss";
   //totals 
   for(std::map<char *, cxxISolutionComp, CHARSTAR_LESS>::iterator iter = this->comps.begin(); iter != this->comps.end(); ++iter)
   {
@@ -930,3 +896,31 @@ void cxxISolution::ORCH_write_output_vars(std::ostream &outstream)
   outstream << std::endl;
   return;
 }
+void cxxISolution::ORCH_write(std::ostream &chemistry_dat, std::ostream &input_dat, std::ostream &output_dat)
+{
+  //
+  // Write orchestra chemistry file definition
+  //
+  chemistry_dat << std::endl << "// Write ORCHESTRA chemistry definitions" << std::endl;
+  // mark for Orchestra include 
+  chemistry_dat << std::endl << "@class: species_reactions () {" << std::endl;
+  this->ORCH_write_chemistry(chemistry_dat);
+  // end orchestra include block
+  chemistry_dat << std::endl << "}" << std::endl;
+  //
+  // Write orchestra input file definition
+  //
+  input_dat << std::endl << "@class: input_file_data () {" << std::endl;
+  this->ORCH_write_input(input_dat);
+  input_dat << std::endl << "}" << std::endl;
+
+  //
+  // Write orchestra output file definition
+  //
+  output_dat << std::endl << "Output_every: 1" << std::endl;
+  this->ORCH_write_output_vars(output_dat);
+
+  //write data to stderr
+  //std::cerr << chemistry_dat.str() << input_dat.str() << output_dat.str();
+
+} 
diff --git a/ISolution.h b/ISolution.h
index b629f013..f747ccd4 100644
--- a/ISolution.h
+++ b/ISolution.h
@@ -41,6 +41,7 @@ public:
         //void dump_xml(std::ostream& os, unsigned int indent = 0)const;
 
 	void ConvertUnits();
+	void ORCH_write(std::ostream &chemistry_dat, std::ostream &input_dat, std::ostream &output_dat);
 	void ORCH_write_chemistry(std::ostream &chemistry_dat);
 	void ORCH_write_input(std::ostream &input_dat);
 	void ORCH_write_output_vars(std::ostream &input_dat);
@@ -61,7 +62,6 @@ private:
   void ORCH_write_chemistry_primary(std::ostream &chemistry_dat);
   void ORCH_write_chemistry_total_O_H(std::ostream &chemistry_dat);
   void ORCH_write_chemistry_alkalinity(std::ostream &chemistry_dat);
-  void ORCH_write_chemistry_species(std::ostream &chemistry_dat);
   void ORCH_write_chemistry_minerals(std::ostream &chemistry_dat);
 };
 
diff --git a/PPassemblage.cxx b/PPassemblage.cxx
index 9dd8bca2..f09e3a1b 100644
--- a/PPassemblage.cxx
+++ b/PPassemblage.cxx
@@ -285,4 +285,76 @@ void cxxPPassemblage::add(const cxxPPassemblage &addee, double extensive)
   //cxxNameDouble eltList;
   this->eltList.add_extensive(addee.eltList, extensive);
 }
+void cxxPPassemblage::ORCH_write_chemistry(std::ostream &chemistry_dat)
+{
+  chemistry_dat << std::endl << "//********* Mineral equilibria" << std::endl;
+//
+// Write minerals
+//
+  for (std::list<cxxPPassemblageComp>::iterator it = this->ppAssemblageComps.begin(); it != this->ppAssemblageComps.end(); ++it) 
+  {
+    struct phase *phase_ptr;
+    int n;
+    phase_ptr = phase_bsearch(it->get_name(), &n, FALSE);
+    assert (phase_ptr != NULL);
+    std::string phase_name(phase_ptr->name);
+    std::replace(phase_name.begin(), phase_name.end(), '(', '[');
+    std::replace(phase_name.begin(), phase_name.end(), ')', ']');
+    chemistry_dat << "@mineral(" << phase_name << ")" << std::endl;
+    chemistry_dat << "@reaction(" << phase_name << ", " << pow(10.0,  -phase_ptr->lk);
+    struct rxn_token *next_token = phase_ptr->rxn_x->token;
+    next_token++;
+    while (next_token->s != NULL || next_token->name != NULL)
+    {
+      chemistry_dat << ", " << next_token->coef;
+      if (next_token->s != NULL)
+      {
+	chemistry_dat << ", " << next_token->s->name;
+      } else {
+	chemistry_dat << ", " << next_token->name;
+      }
+      next_token++;
+    }
+    chemistry_dat << ")" << std::endl;
+  }
+}
 
+void cxxPPassemblage::ORCH_write_output_vars(std::ostream &outstream)
+{
+  //
+  //  Write ppassemblage output variables
+  //
+  outstream << "\tstart_ppassemblage";
+  for (std::list<cxxPPassemblageComp>::iterator it = this->ppAssemblageComps.begin(); it != this->ppAssemblageComps.end(); ++it) 
+  {
+    outstream << "\t" << it->get_name() << ".min";
+  }
+  outstream << "\tend_ppassemblage";
+}
+void cxxPPassemblage::ORCH_read(std::vector <std::pair <std::string, double> > output_vector, std::vector < std::pair < std::string, double > >::iterator &it)
+{
+  while (it->first.compare("end_ppassemblage") != 0) 
+  {
+    it++;
+  } 
+}
+
+void cxxPPassemblage::ORCH_store_global(std::map < std::string, double > output_map)
+{
+  int i;
+  std::map < std::string, double >::iterator it;
+  for (i = 0; i < count_unknowns; i++)
+  {
+    // MB, ALK, CB, SOLUTION_PHASE_BOUNDARY, MU, AH2O
+    switch (x[i]->type)
+    {
+      case PP:
+	std::string name(x[i]->phase->name);
+	name.append(".min");
+	it = output_map.find(name);
+	assert (it != output_map.end());
+	x[i]->moles = it->second;
+	break;
+    }
+  }
+}
diff --git a/PPassemblage.h b/PPassemblage.h
index 659e8304..d379bada 100644
--- a/PPassemblage.h
+++ b/PPassemblage.h
@@ -31,23 +31,27 @@ public:
 
         void read_raw(CParser& parser);
 
+	const cxxNameDouble& get_totals()const
+	{
+	  return this->totals;
+	};
+
 #ifdef USE_MPI
 	void mpi_pack(std::vector<int>& ints, std::vector<double>& doubles);
 	void mpi_unpack(int *ints, int *ii, double *doubles, int *dd);
 #endif
 
 	void totalize();
+	void ORCH_write_chemistry(std::ostream &chemistry_dat);
+	void ORCH_write_output_vars(std::ostream &outstream);
+	void ORCH_read(std::vector <std::pair <std::string, double> > output_vector, std::vector < std::pair < std::string, double > >::iterator &it);
+	void ORCH_store_global(std::map < std::string, double > output_map);
 
 private:
 	void add(const cxxPPassemblage &addee, double extensive);
 	// not written
         void dump_xml(std::ostream& os, unsigned int indent = 0)const;
 
-	const cxxNameDouble& get_totals()const
-	{
-	  return this->totals;
-	};
-
 protected:
         std::list<cxxPPassemblageComp> ppAssemblageComps;
         cxxNameDouble eltList;
diff --git a/Solution.cxx b/Solution.cxx
index 04cc70e5..b3b866f8 100644
--- a/Solution.cxx
+++ b/Solution.cxx
@@ -976,9 +976,11 @@ void cxxSolution::ORCH_read(std::vector <std::pair <std::string, double> > outpu
   this->total_o = it->second; it++;
   this->cb = it->second; it++;
   this->mass_water = it->second * gfw_water; it++;
+  this->mass_water = 1.0;
   this->total_alkalinity = it->second; it++;
   it++; //orch total H+
   it++; //orch total e-
+  it++; //orch total H2O
   //cxxNameDouble totals;
   char token[MAX_LENGTH];
   while (it->first.compare("end_totals") != 0) 
diff --git a/Solution.h b/Solution.h
index ef519e4d..a0efe286 100644
--- a/Solution.h
+++ b/Solution.h
@@ -62,10 +62,11 @@ public:
 	double get_total_element(char *string)const;
 	void set_total(char *string, double value);
 
+	const cxxNameDouble& get_totals(void)const {return this->totals;}
 	void set_totals(cxxNameDouble &nd) {this->totals = nd; this->totals.type = cxxNameDouble::ND_ELT_MOLES;}
 	void clear_totals() {this->totals.clear();}
 
-
+	const cxxNameDouble& get_master_activity(void)const {return this->master_activity;}
 	double get_master_activity(char *string)const;
 	void set_master_activity(char *string, double value);
 
diff --git a/StorageBin.cxx b/StorageBin.cxx
index 4b8b3064..71bee254 100644
--- a/StorageBin.cxx
+++ b/StorageBin.cxx
@@ -9,7 +9,7 @@
 #include <iostream>     // std::cout std::cerr
 #include "Utils.h"   // define first
 #include "StorageBin.h"
-#include "Solution.h"
+#include "System.h"
 #define EXTERNAL extern
 #include "global.h"
 #include "phqalloc.h"
@@ -25,7 +25,73 @@ cxxStorageBin::cxxStorageBin()
 {
         // default constructor for cxxStorageBin 
 }
+cxxStorageBin::cxxStorageBin(struct Use *use_ptr)
+{
+  //Construct from use pointer
 
+  // Solution
+  if (use_ptr->solution_ptr != NULL)
+  {
+    cxxSolution entity(use_ptr->solution_ptr);
+    this->setSolution(use_ptr->n_solution_user, &entity);
+  }
+  // Exchange
+  if (use_ptr->exchange_ptr != NULL)
+  {
+    cxxExchange entity(use_ptr->exchange_ptr);
+    this->setExchange(use_ptr->n_exchange_user, &entity);
+  }	
+  // gas_phase
+  if (use_ptr->gas_phase_ptr != NULL)
+  {
+    cxxGasPhase entity(use_ptr->gas_phase_ptr);
+    this->setGasPhase(use_ptr->n_gas_phase_user, &entity);
+  }
+  // kinetics
+  if (use_ptr->kinetics_ptr != NULL)
+  {
+    cxxKinetics entity(use_ptr->kinetics_ptr);
+    this->setKinetics(use_ptr->n_kinetics_user, &entity);
+  }
+  // pp_assemblage
+  if (use_ptr->pp_assemblage_ptr != NULL)
+  {
+    cxxPPassemblage entity(use_ptr->pp_assemblage_ptr);
+    this->setPPassemblage(use_ptr->n_pp_assemblage_user, &entity);
+  }	
+  // s_s_assemblage
+  if (use_ptr->s_s_assemblage_ptr != NULL)
+  {
+    cxxSSassemblage entity(use_ptr->s_s_assemblage_ptr);
+    this->setSSassemblage(use_ptr->n_s_s_assemblage_user, &entity);
+  }	
+  // surface
+  if (use_ptr->surface_ptr != NULL)
+  {
+    cxxSurface entity(use_ptr->surface_ptr);
+    this->setSurface(use_ptr->n_surface_user, &entity);
+  }	
+  // mix
+  if (use_ptr->mix_ptr != NULL)
+  {
+    cxxMix entity(use_ptr->mix_ptr);
+    this->setMix(use_ptr->n_mix_user, &entity);
+  }	
+  // reaction
+  if (use_ptr->irrev_ptr != NULL)
+  {
+    cxxReaction entity(use_ptr->irrev_ptr);
+    this->setReaction(use_ptr->n_irrev_user, &entity);
+  }	
+  // reaction temperature
+  if (use_ptr->temperature_ptr != NULL)
+  {
+    cxxTemperature entity(use_ptr->temperature_ptr);
+    this->setTemperature(use_ptr->n_temperature_user, &entity);
+  }
+  // set system
+  this->setSystem(use_ptr);
+}
 cxxStorageBin::~cxxStorageBin()
 {
 }
@@ -917,3 +983,107 @@ cxxExchange *cxxStorageBin::mix_cxxExchange(cxxMix &mixmap)
 	return(new_exch_ptr);
 }
 #endif
+
+void cxxStorageBin::ORCH_write(std::ostream &chemistry_dat, std::ostream &input_dat, std::ostream &output_dat)
+{
+  //std::ostringstream oss;
+
+  // Liter concentrations
+  this->system.ORCH_write(chemistry_dat, input_dat, output_dat);
+
+}
+void cxxStorageBin::setSystem(struct Use *use_ptr)
+{
+  // Initialize
+  this->system.Initialize();
+  // Solution
+  if (use_ptr->solution_ptr != NULL)
+  {
+    std::map <int, cxxSolution>::iterator it = this->Solutions.find(use_ptr->n_solution_user);
+    if (it != this->Solutions.end())
+    {
+      this->system.setSolution(&(it->second));
+    }
+  }
+  // Exchange
+  if (use_ptr->exchange_ptr != NULL)
+  {
+    std::map <int, cxxExchange>::iterator it = this->Exchangers.find(use_ptr->n_exchange_user);
+    if (it != this->Exchangers.end())
+    {
+      this->system.setExchange(&(it->second));
+    }
+   }	
+  // gas_phase
+  if (use_ptr->gas_phase_ptr != NULL)
+  {
+    std::map <int, cxxGasPhase>::iterator it = this->GasPhases.find(use_ptr->n_gas_phase_user);
+    if (it != this->GasPhases.end())
+    {
+      this->system.setGasPhase(&(it->second));
+    }
+  }
+  // kinetics
+  if (use_ptr->kinetics_ptr != NULL)
+  {
+    std::map <int, cxxKinetics>::iterator it = this->Kinetics.find(use_ptr->n_kinetics_user);
+    if (it != this->Kinetics.end())
+    {
+      this->system.setKinetics(&(it->second));
+    }
+  }
+  // pp_assemblage
+  if (use_ptr->pp_assemblage_ptr != NULL)
+  {
+    std::map <int, cxxPPassemblage>::iterator it = this->PPassemblages.find(use_ptr->n_pp_assemblage_user);
+    if (it != this->PPassemblages.end())
+    {
+      this->system.setPPassemblage(&(it->second));
+    }
+  }	
+  // s_s_assemblage
+  if (use_ptr->s_s_assemblage_ptr != NULL)
+  {
+    std::map <int, cxxSSassemblage>::iterator it = this->SSassemblages.find(use_ptr->n_s_s_assemblage_user);
+    if (it != this->SSassemblages.end())
+    {
+      this->system.setSSassemblage(&(it->second));
+    }
+  }	
+  // surface
+  if (use_ptr->surface_ptr != NULL)
+  {
+    std::map <int, cxxSurface>::iterator it = this->Surfaces.find(use_ptr->n_surface_user);
+    if (it != this->Surfaces.end())
+    {
+      this->system.setSurface(&(it->second));
+    }
+  }	
+  // mix
+  if (use_ptr->mix_ptr != NULL)
+  {
+    std::map <int, cxxMix>::iterator it = this->Mixes.find(use_ptr->n_mix_user);
+    if (it != this->Mixes.end())
+    {
+      this->system.setMix(&(it->second));
+    }
+  }	
+  // reaction
+  if (use_ptr->irrev_ptr != NULL)
+  {
+    std::map <int, cxxReaction>::iterator it = this->Reactions.find(use_ptr->n_irrev_user);
+    if (it != this->Reactions.end())
+    {
+      this->system.setReaction(&(it->second));
+    }
+  }	
+  // reaction temperature
+  if (use_ptr->temperature_ptr != NULL)
+  {
+    std::map <int, cxxTemperature>::iterator it = this->Temperatures.find(use_ptr->n_temperature_user);
+    if (it != this->Temperatures.end())
+    {
+      this->system.setTemperature(&(it->second));
+    }
+  }
+}
\ No newline at end of file
diff --git a/StorageBin.h b/StorageBin.h
index a3ee4123..22d568cc 100644
--- a/StorageBin.h
+++ b/StorageBin.h
@@ -10,6 +10,7 @@
 #include "PPassemblage.h"
 #include "SSassemblage.h"
 #include "Surface.h"
+#include "System.h"
 #include "cxxMix.h"
 #include "Reaction.h"
 #include "Temperature.h"
@@ -25,7 +26,7 @@ class cxxStorageBin
 
 public:
         cxxStorageBin();
-
+        cxxStorageBin(struct Use *use_ptr);
         ~cxxStorageBin();
 
         void import_phreeqc(void);
@@ -136,7 +137,51 @@ public:
 		Surfaces.erase(n_user);
 	}
 
-        void dump_raw(std::ostream& s_oss, unsigned int indent)const;
+	cxxMix *getMix(int n_user) { 
+		if (this->Mixes.find(n_user) != this->Mixes.end()) {
+			return(&(this->Mixes.find(n_user)->second));
+		} 
+		return (NULL);
+	}
+	void setMix(int n_user, cxxMix *entity) {
+	  if (entity == NULL) return;
+	  Mixes[n_user] = *entity;
+	}
+	void removeMix(int n_user) { 
+		Mixes.erase(n_user);
+	}
+
+ 	cxxReaction *getReaction(int n_user) { 
+		if (this->Reactions.find(n_user) != this->Reactions.end()) {
+			return(&(this->Reactions.find(n_user)->second));
+		} 
+		return (NULL);
+	}
+	void setReaction(int n_user, cxxReaction *entity) {
+	  if (entity == NULL) return;
+	  Reactions[n_user] = *entity;
+	}
+	void removeReaction(int n_user) { 
+		Reactions.erase(n_user);
+	}
+
+ 	cxxTemperature *getTemperature(int n_user) { 
+		if (this->Temperatures.find(n_user) != this->Temperatures.end()) {
+			return(&(this->Temperatures.find(n_user)->second));
+		} 
+		return (NULL);
+	}
+	void setTemperature(int n_user, cxxTemperature *entity) {
+	  if (entity == NULL) return;
+	  Temperatures[n_user] = *entity;
+	}
+	void removeTemperature(int n_user) { 
+		Temperatures.erase(n_user);
+	}
+
+	void setSystem(struct Use *use_ptr); 
+
+	void dump_raw(std::ostream& s_oss, unsigned int indent)const;
 
         void dump_raw(std::ostream& s_oss, int i, unsigned int indent);
 
@@ -161,7 +206,7 @@ public:
 	void mpi_send(int n, int task_number);
 	void mpi_recv(int task_number);
 #endif
-
+	void ORCH_write(std::ostream &chemistry_dat, std::ostream &input_dat, std::ostream &output_dat);
 protected:
 	// Tidied classes
 	std::map<int, cxxSolution>      Solutions;
@@ -180,7 +225,7 @@ protected:
 	std::map<int, cxxMix>           Mixes;
 	std::map<int, cxxReaction>      Reactions;
 	std::map<int, cxxTemperature>   Temperatures;
-
+	cxxSystem                       system;
 public:
         //static std::map<int, cxxStorageBin>& map;
 
diff --git a/Surface.cxx b/Surface.cxx
index a859cd40..81c6565d 100644
--- a/Surface.cxx
+++ b/Surface.cxx
@@ -680,6 +680,7 @@ void cxxSurface::totalize()
   for (std::list<cxxSurfaceComp>::const_iterator it = surfaceComps.begin(); it != surfaceComps.end(); ++it) 
   {
     this->totals.add_extensive(it->get_totals(), 1.0);
+    this->totals.add("Charge", it->get_charge_balance());
   }
   return;
 }
diff --git a/SurfaceComp.h b/SurfaceComp.h
index 4eec084f..67a4cce0 100644
--- a/SurfaceComp.h
+++ b/SurfaceComp.h
@@ -27,6 +27,8 @@ public:
         char *get_phase_name()const {return this->phase_name;}
         char *get_rate_name()const {return this->rate_name;}
         char *get_formula()const {return this->formula;}
+        double get_charge_balance()const        { return this->charge_balance;}
+        void   set_charge_balance(double d)     { this->charge_balance = d;}
 
         static struct surface_comp *cxxSurfaceComp2surface_comp(std::list<cxxSurfaceComp>& el);
 
diff --git a/System.cxx b/System.cxx
new file mode 100644
index 00000000..ed5e3c54
--- /dev/null
+++ b/System.cxx
@@ -0,0 +1,513 @@
+#include "System.h"
+#include <algorithm>   // std::replace 
+extern void ORCH_write_chemistry_species(std::ostream &chemistry_dat);
+cxxSystem::cxxSystem(void)
+{
+  this->solution = NULL;
+  this->exchange = NULL;
+  this->ppassemblage = NULL;
+  this->gasphase = NULL;
+  this->ssassemblage = NULL;
+  this->kinetics = NULL;
+  this->surface = NULL;
+  this->mix = NULL;
+  this->reaction = NULL;
+  this->temperature = NULL;
+}
+
+cxxSystem::~cxxSystem(void)
+{
+}
+
+void cxxSystem::Initialize(void)
+{
+  this->solution = NULL;
+  this->exchange = NULL;
+  this->ppassemblage = NULL;
+  this->gasphase = NULL;
+  this->ssassemblage = NULL;
+  this->kinetics = NULL;
+  this->surface = NULL;
+  this->mix = NULL;
+  this->reaction = NULL;
+  this->temperature = NULL;
+}
+void cxxSystem::totalize(void)
+{
+  //initialize
+  this->totals.clear();
+  //add solution
+  if (this->solution != NULL)
+  {
+    this->totals["O"] = this->solution->get_total_o();
+    this->totals["H"] = this->solution->get_total_h();
+    this->totals["Charge"] = this->solution->get_cb();
+    this->totals.add_extensive(this->solution->get_totals(), 1.0);
+  }
+  if (this->exchange != NULL)
+  {
+    this->exchange->totalize();
+    this->totals.add_extensive(this->exchange->get_totals(), 1.0);
+  }
+  if (this->ppassemblage != NULL)
+  {
+    this->ppassemblage->totalize();
+    this->totals.add_extensive(this->ppassemblage->get_totals(), 1.0);
+  }
+  if (this->gasphase != NULL)
+  {
+    this->gasphase->totalize();
+    this->totals.add_extensive(this->gasphase->get_totals(), 1.0);
+  }
+  if (this->ssassemblage != NULL)
+  {
+    this->ssassemblage->totalize();
+    this->totals.add_extensive(this->ssassemblage->get_totals(), 1.0);
+  }
+  if (this->surface != NULL)
+  {
+    this->ssassemblage->totalize();
+    this->totals.add_extensive(this->surface->get_totals(), 1.0);
+  }
+
+  //Need to handle the following 3 reactions:
+  //
+  //this->kinetics = NULL;
+  //this->mix = NULL;
+  //this->reaction = NULL;
+  //this->temperature = NULL;
+  //this->totals.dump_raw(std::cerr, 1);
+  return;
+}
+void cxxSystem::ORCH_write(std::ostream &chemistry_dat, std::ostream &input_dat, std::ostream &output_dat)
+{
+  //
+  // chemistry file
+  //
+  this->totalize();
+  // calculate orchestra components from PHREEQC totals
+  this->ORCH_components();
+  // write definitions for chemistry
+  chemistry_dat << std::endl << "@class: species_reactions () {" << std::endl;
+  this->ORCH_write_chemistry_water(chemistry_dat);
+  this->ORCH_write_chemistry_primary(chemistry_dat);
+  this->ORCH_write_chemistry_total_O_H(chemistry_dat);
+  ORCH_write_chemistry_species(chemistry_dat);
+
+  // add definitions for pure phases
+  if (this->ppassemblage != NULL) this->ppassemblage->ORCH_write_chemistry(chemistry_dat);
+
+  // add definitions for other PHREEQC entities here
+
+  // finish up
+  chemistry_dat << std::endl << "}" << std::endl;
+  //
+  // input file
+  //
+  input_dat << std::endl << "@class: input_file_data () {" << std::endl;
+  this->ORCH_write_input(input_dat);
+  input_dat << std::endl << "}" << std::endl;
+  //
+  // output file
+  //
+  output_dat << std::endl << "Output_every: 1" << std::endl;
+  output_dat << "Var:";
+  this->ORCH_write_output_vars(output_dat);
+  // add definitions for pure phases
+  if (this->ppassemblage != NULL) this->ppassemblage->ORCH_write_output_vars(output_dat);
+
+  //finish up
+  output_dat << std::endl;
+};
+
+void cxxSystem::ORCH_write_chemistry_water(std::ostream &chemistry_dat)
+{
+  //
+  //  Write water entities
+  //
+  chemistry_dat << std::endl << "//********* The water entities" << std::endl;
+
+  //  e-, total hydrogen
+  chemistry_dat << "@entity(e-, diss, 0)" << std::endl;
+  chemistry_dat << "@Calc: (1, \"e-.act = 10.^(-pe)\")" << std::endl;
+  chemistry_dat << "@solve (pe, 1e-6, lin, 1, e-.liter,  1e-14)" << std::endl;
+
+  //  H+, charge balance
+  chemistry_dat << "@Calc: (1, \"H+.act = 10.^(-pH)\")" << std::endl;
+  chemistry_dat << "@solve (pH, 1e-6, lin, 1, H+.liter,  1e-14)" << std::endl;
+ 
+  // H2O
+  chemistry_dat << "@entity(" << s_h2o->name << ", diss, 55.506)" << std::endl;
+  chemistry_dat << std::endl;
+}
+
+void cxxSystem::ORCH_write_chemistry_primary(std::ostream &chemistry_dat)
+{
+  chemistry_dat << std::endl << "//********* The primary species" << std::endl;
+  //
+  //  Write other master species definitions, i.e. primary entities
+  //
+  chemistry_dat << "@species(H+, 1)" << std::endl;
+  for(cxxNameDouble::iterator iter = this->totals.begin() ; iter != this->totals.end(); ++iter)
+  {
+    std::string name(iter->first);
+    if (name == "H(1)" || name == "E" || name == "H" || name == "O" || name == "Charge") continue;
+    struct element *elt;
+    char *element_name = string_hsave(name.c_str());
+    elt = element_store(element_name);
+    assert(elt != NULL);
+    struct species *s_ptr;
+    s_ptr = elt->master->s;
+    assert(s_ptr != NULL);
+    chemistry_dat << "@species(" << s_ptr->name << ", " << s_ptr->z << ")" << std::endl;
+    // regular mole balance
+    chemistry_dat << "@primary_entity(" << s_ptr->name << ", 1e-9, liter, 1.0e-9)" << std::endl;
+  } 
+  chemistry_dat << std::endl;
+}
+
+void cxxSystem::ORCH_write_chemistry_total_O_H(std::ostream &chemistry_dat)
+{
+  chemistry_dat << std::endl << "//********* Define total hydrogen and oxygen" << std::endl;
+
+  // Write H total equation
+  chemistry_dat << "@var: total_hydrogen 0" << std::endl;
+  chemistry_dat << "@calc: (5, \"total_hydrogen = 0";
+  int i;
+  for (i = 0; i < count_s_x; i++) 
+  {
+    // write in terms of orchestra components
+    if (s_x[i]->primary != NULL || (s_x[i]->secondary != NULL && s_x[i]->secondary->in == TRUE))
+    {
+      if (s_x[i]->h != 0) 
+      {
+	chemistry_dat << "+";
+	if (s_x[i]->h != 1)
+	{
+	  chemistry_dat << s_x[i]->h << "*";
+	}
+	chemistry_dat << "{" << s_x[i]->name << ".liter}";
+      }
+    }
+  } 
+  chemistry_dat << "\")" << std::endl;
+
+  // Write O total equation
+  chemistry_dat << "@var: total_oxygen 0" << std::endl;
+  chemistry_dat << "@calc: (5, \"total_oxygen = 0";
+  for (i = 0; i < count_s_x; i++) 
+  {
+    if (s_x[i]->o != 0) 
+    {
+      // write in terms of orchestra components
+      if (s_x[i]->primary != NULL || (s_x[i]->secondary != NULL && s_x[i]->secondary->in == TRUE))
+      {
+	chemistry_dat << "+";
+	if (s_x[i]->o != 1)
+	{
+	  chemistry_dat << s_x[i]->o << "*";
+	}
+	chemistry_dat << "{" << s_x[i]->name << ".liter}";
+      }
+    }
+  } 
+  chemistry_dat << "\")" << std::endl;
+  chemistry_dat << std::endl;
+}
+#ifdef SKIP
+void cxxSystem::ORCH_write_input(std::ostream &input_dat)
+{
+
+
+//
+//  Write orchestra input file info
+//
+  std::ostringstream headings, data;
+  data.precision(DBL_DIG - 1);
+  headings << "var:   ";
+  data << "data:  ";
+
+
+  headings << "tempc\t";
+  data   << this->solution->get_tc() << "\t";
+
+  headings << "pH\t";
+  data   << this->solution->get_ph() << "\t";
+
+  headings << "pe\t";
+  data   << this->solution->get_pe() << "\t";
+
+  headings << "H2O.act\t";
+  data   << 1 << "\t";
+
+  for (cxxNameDouble::iterator iter = this->totals.begin(); iter != this->totals.end(); ++iter) 
+  {
+    if (iter->first == "O") continue;
+    std::string master_name;
+    struct master *master_ptr;
+    master_ptr = master_bsearch (iter->first);
+    assert (master_ptr != NULL);
+    double coef = master_ptr->coef;
+    if (master_ptr->coef == 0) 
+    {
+      coef = 1;
+    }
+    if (iter->first == "H")
+    {
+      headings << "system_hydrogen"  << "\t";
+      data << iter->second << "\t";
+      continue;
+    } else
+    {
+      headings << master_ptr->s->name << ".liter"  << "\t";
+      data << iter->second / coef << "\t";
+    }
+
+    //  activity estimate
+    cxxNameDouble ma = this->solution->get_master_activity();
+    cxxNameDouble::iterator it = ma.find(iter->first);
+    if (it == ma.end())
+    {
+      it = ma.find(master_ptr->s->secondary->elt->name);
+    }
+    headings << master_ptr->s->name << ".act\t";
+    if (it != ma.end())
+    {
+      data << pow(10., it->second) << "\t";
+    } else
+    {
+      data << 1e-9 << "\t";
+    }
+  }
+  // Isotopes
+  //s_oss << "-Isotopes" << std::endl;
+  /*
+  {
+          for (std::list<cxxSolutionIsotope>::const_iterator it = this->isotopes.begin(); it != isotopes.end(); ++it) {
+                  it->dump_raw(s_oss, indent + 2);
+          }
+  }
+  */
+
+  // Write data to string
+  input_dat << headings.str() << std::endl;
+  input_dat << data.str() << std::endl;
+
+  return;
+}
+#endif
+void cxxSystem::ORCH_write_input(std::ostream &input_dat)
+{
+
+
+//
+//  Write orchestra input file info
+//
+  std::ostringstream headings, data;
+  data.precision(DBL_DIG - 1);
+  headings << "var:   ";
+  data << "data:  ";
+
+  headings << "tempc\t";
+  data   << this->solution->get_tc() << "\t";
+
+  headings << "pH\t";
+  data   << this->solution->get_ph() << "\t";
+
+  headings << "pe\t";
+  data   << this->solution->get_pe() << "\t";
+
+  headings << "H2O.act\t";
+  data   << 1 << "\t";
+
+  headings << "I\t";
+  data   << this->solution->get_mu() << "\t";
+
+  for (cxxNameDouble::iterator iter = this->orch_totals.begin(); iter != this->orch_totals.end(); ++iter) 
+  {
+    headings << iter->first << ".liter"  << "\t";
+    data << iter->second << "\t";
+  }
+
+  //  activity estimate
+  for (cxxNameDouble::iterator iter = this->totals.begin(); iter != this->totals.end(); ++iter) 
+  {
+    if (iter->first == "O" || iter->first == "Charge") continue;
+    std::string master_name;
+    struct master *master_ptr;
+    master_ptr = master_bsearch (iter->first);
+    assert (master_ptr != NULL);
+  
+    cxxNameDouble ma = this->solution->get_master_activity();
+    cxxNameDouble::iterator it = ma.find(iter->first);
+    if (it == ma.end())
+    {
+      it = ma.find(master_ptr->s->secondary->elt->name);
+    }
+    headings << master_ptr->s->name << ".act\t";
+    if (it != ma.end())
+    {
+      data << pow(10., it->second) << "\t";
+    } else
+    {
+      data << 1e-9 << "\t";
+    }
+  }
+
+  // Isotopes
+  //s_oss << "-Isotopes" << std::endl;
+  /*
+  {
+          for (std::list<cxxSolutionIsotope>::const_iterator it = this->isotopes.begin(); it != isotopes.end(); ++it) {
+                  it->dump_raw(s_oss, indent + 2);
+          }
+  }
+  */
+
+  // Write data to string
+  input_dat << headings.str() << std::endl;
+  input_dat << data.str() << std::endl;
+
+  return;
+}
+
+void cxxSystem::ORCH_components(void)
+{
+  // translate from H, O, charge to H+tot, e-tot, and H2Otot
+  cxxNameDouble::iterator it;
+  cxxNameDouble temp_totals;
+  //
+  // Set names in orch_totals to master species names
+  //
+  this->orch_totals.clear();
+  temp_totals = this->totals;
+  temp_totals.erase("H");
+  temp_totals.erase("O");
+  temp_totals.erase("Charge");
+  //temp_totals.dump_raw(std::cerr, 1);
+  for (it = temp_totals.begin(); it != temp_totals.end(); it++) 
+  {
+    struct element *elt_ptr; 
+    elt_ptr = element_store(it->first);
+    assert(elt_ptr != NULL);
+    struct master *master_ptr;
+    master_ptr = elt_ptr->primary;
+    assert(master_ptr != NULL);
+    double coef = master_ptr->coef;
+    if (master_ptr->coef == 0) 
+    {
+      coef = 1;
+    }
+    this->orch_totals[master_ptr->s->name] = it->second / coef;
+  }
+  //
+  // Calculate h2otot
+  //
+  it = this->totals.find("O");
+  assert (it != this->totals.end());
+  double h2otot = it->second;
+  // subtract O in master species
+  for (it = temp_totals.begin(); it != temp_totals.end(); it++)
+  {
+    struct element *elt_ptr; 
+    elt_ptr = element_store(it->first);
+    struct master *master_ptr;
+    master_ptr = elt_ptr->primary;
+    h2otot -= it->second * master_ptr->s->o;
+  }
+  //
+  // Calculate htot
+  //
+  it = this->totals.find("H");
+  assert (it != this->totals.end());
+  double htot = it->second - 2*h2otot;
+  // subtract O in master species
+  for (it = temp_totals.begin(); it != temp_totals.end(); it++)
+  {
+    struct element *elt_ptr; 
+    elt_ptr = element_store(it->first);
+    struct master *master_ptr;
+    master_ptr = elt_ptr->primary;
+    htot -= it->second * master_ptr->s->h;
+  }
+  //
+  // Calculate etot
+  //
+  it = this->totals.find("Charge");
+  assert (it != this->totals.end());
+  double etot = -it->second + htot;
+  // subtract O in master species
+  for (it = temp_totals.begin(); it != temp_totals.end(); it++)
+  {
+    struct element *elt_ptr; 
+    elt_ptr = element_store(it->first);
+    struct master *master_ptr;
+    master_ptr = elt_ptr->primary;
+    etot += it->second * master_ptr->s->z;
+  }
+  //
+  // store h2otot, htot, etot in orch_totals
+  //
+  this->orch_totals["H2O"] = h2otot;
+  this->orch_totals["H+"] = htot;
+  this->orch_totals["e-"] = etot;
+  this->orch_totals.dump_raw(std::cerr, 1);
+}
+void cxxSystem::ORCH_write_output_vars(std::ostream &outstream)
+{
+  outstream << "\tnr_iter";
+  //
+  //  Serialize solution
+  //
+  outstream << "\tstart_solution";
+  //tc
+  outstream << "\ttempc";
+  //pH
+  outstream << "\tpH";
+  //pe
+  outstream << "\tpe";
+  //mu
+  outstream << "\tI";
+  //ah2o
+  outstream << "\tH2O.act";
+  //total_h;
+  outstream << "\ttotal_hydrogen";
+  //total_o;
+  outstream << "\ttotal_oxygen";
+  //cb
+  outstream << "\tchargebalance";
+  //mass_water;
+  outstream << "\tH2O.con";
+  //total_alkalinity;
+  outstream << "\tAlkalinity";
+  //orchestra master variables
+  outstream << "\tH+.diss";
+  outstream << "\te-.diss";
+  outstream << "\tH2O.diss";
+  //
+  // Write totals
+  for (cxxNameDouble::iterator it = this->orch_totals.begin(); it != this->orch_totals.end(); it++)
+  {
+    if ( it->first == "H+" || it->first == "H2O" || it->first == "e-") continue;
+    outstream << "\t" << it->first << ".diss";
+  }
+  outstream << "\tend_totals";
+  for (cxxNameDouble::iterator it = this->orch_totals.begin(); it != this->orch_totals.end(); it++)
+  {
+    if ( it->first == "H+" || it->first == "H2O" || it->first == "e-") continue;
+    struct species *s_ptr = s_search(it->first);
+    outstream << "\t" << it->first << ".act";
+  }
+  outstream << "\tend_master_activities";
+  //
+  // Write all species activities and concentrations
+  //
+  int i;
+  for (i = 0; i < count_s_x; i++) {
+    std::string name(s_x[i]->name);
+    std::replace(name.begin(), name.end(), '(', '[');
+    std::replace(name.begin(), name.end(), ')', ']');
+    outstream << "\t" << name.c_str() << ".act" << "\t" << name.c_str() << ".con";
+  }
+  outstream << "\tend_species";
+}
\ No newline at end of file
diff --git a/System.h b/System.h
new file mode 100644
index 00000000..6dfe7bb9
--- /dev/null
+++ b/System.h
@@ -0,0 +1,73 @@
+#if !defined(SYSTEM_H_INCLUDED)
+#define SYSTEM_H_INCLUDED
+#include "Solution.h"
+#include "Exchange.h"
+#include "GasPhase.h"
+#include "cxxKinetics.h"
+#include "PPassemblage.h"
+#include "SSassemblage.h"
+#include "Surface.h"
+#include "cxxMix.h"
+#include "Reaction.h"
+#include "Temperature.h"
+class cxxSystem
+{
+public:
+  cxxSystem(void);
+public:
+  ~cxxSystem(void);
+  void Initialize(void);
+  void setSolution(cxxSolution *entity) { 
+      this->solution = entity;
+  }
+  void setExchange(cxxExchange *entity) { 
+    this->exchange = entity;
+  }
+  void setPPassemblage(cxxPPassemblage *entity) { 
+    this->ppassemblage = entity;
+  }
+  void setGasPhase(cxxGasPhase *entity) { 
+    this->gasphase = entity;
+  }
+  void setSSassemblage(cxxSSassemblage *entity) { 
+    this->ssassemblage = entity;
+  }
+  void setKinetics(cxxKinetics *entity) { 
+    this->kinetics = entity;
+  }
+  void setSurface(cxxSurface *entity) { 
+    this->surface = entity;
+  }
+  void setMix(cxxMix *entity) { 
+    this->mix = entity;
+  }
+  void setReaction(cxxReaction *entity) { 
+    this->reaction = entity;
+  }
+  void setTemperature(cxxTemperature *entity) { 
+    this->temperature = entity;
+  }
+  void totalize();
+  void ORCH_components();
+  void ORCH_write(std::ostream &chemistry_dat, std::ostream &input_dat, std::ostream &output_dat);
+  void ORCH_write_chemistry_water(std::ostream &chemistry_dat);
+  void ORCH_write_chemistry_primary(std::ostream &chemistry_dat);
+  void ORCH_write_chemistry_total_O_H(std::ostream &chemistry_dat);
+  void ORCH_write_output_vars(std::ostream &outstream);
+  void ORCH_write_input(std::ostream &input_dat);
+
+private:
+  cxxSolution *solution;
+  cxxExchange *exchange;
+  cxxPPassemblage *ppassemblage;
+  cxxGasPhase *gasphase;
+  cxxSSassemblage *ssassemblage;
+  cxxKinetics *kinetics;
+  cxxSurface *surface;
+  cxxMix *mix;
+  cxxReaction *reaction;
+  cxxTemperature *temperature;
+  cxxNameDouble totals;
+  cxxNameDouble orch_totals;
+};
+#endif // !defined(SYSTEM_H_INCLUDED)