First stab at aqueous kinetics. This has not been checked.

Cantera/src/kinetics/AqueousKinetics.cpp

@@ -0,0 +1,550 @@
+/**
+ *  @file AqueousKinetics.cpp 
+ *
+ * Homogeneous kinetics in an aqueous phase, either condensed
+ * or dilute in salts
+ * 
+ */
+/*
+ * Copywrite (2006) Sandia Corporation. Under the terms of
+ * Contract DE-AC04-94AL85000 with Sandia Corporation, the
+ * U.S. Government retains certain rights in this software.
+ */
+/*
+ *  $Date$
+ *  $Revision$
+ */
+
+
+
+#ifdef WIN32
+#pragma warning(disable:4786)
+#pragma warning(disable:4503)
+#endif
+
+#include "AqueousKinetics.h"
+#include "ReactionData.h"
+#include "RateCoeffMgr.h"
+
+#include <iostream>
+using namespace std;
+
+namespace Cantera {
+
+  /**
+   * Construct an empty reaction mechanism.
+   */    
+  AqueousKinetics::AqueousKinetics(thermo_t* thermo) :
+    Kinetics(),
+    m_kk(0), 
+    m_nfall(0), 
+    m_nirrev(0), 
+    m_nrev(0),
+    m_finalized(false)
+  {
+    if (thermo != 0) addPhase(*thermo);
+    m_kdata = new AqueousKineticsData;
+    m_kdata->m_temp = 0.0;
+    m_rxnstoich = new ReactionStoichMgr;
+  }
+
+  AqueousKinetics::~AqueousKinetics() {
+    delete m_kdata; 
+    delete m_rxnstoich;
+  }
+
+  /**
+   * Update temperature-dependent portions of reaction rates and
+   * falloff functions.
+   */
+  void AqueousKinetics::
+  update_T() {}
+
+  void AqueousKinetics::
+  update_C() {}
+
+  void AqueousKinetics::_update_rates_T() {
+    doublereal T = thermo().temperature();
+    // m_kdata->m_logStandConc = log(thermo().standardConcentration()); 
+    doublereal logT = log(T);
+    m_rates.update(T, logT, &m_kdata->m_rfn[0]);
+ 
+
+    m_kdata->m_temp = T;
+    updateKc();
+    m_kdata->m_ROP_ok = false;
+   
+  };
+
+
+  /**
+   * Update properties that depend on concentrations. Currently only
+   * the enhanced collision partner concentrations are updated here.
+   */         
+  void AqueousKinetics::
+  _update_rates_C() {
+    thermo().getActivityConcentrations(&m_conc[0]);
+   
+    m_kdata->m_ROP_ok = false;
+  }
+
+  /**
+   * Update the equilibrium constants in molar units.
+   */
+  void AqueousKinetics::updateKc() {
+    int i, irxn;
+    vector_fp& m_rkc = m_kdata->m_rkcn;
+    doublereal rt = GasConstant*   m_kdata->m_temp;
+        
+    thermo().getStandardChemPotentials(&m_grt[0]);
+    fill(m_rkc.begin(), m_rkc.end(), 0.0);
+    int nsp = thermo().nSpecies();
+    for (int k = 0; k < nsp; k++) {
+      doublereal logStandConc_k = thermo().logStandardConc(k);
+      m_grt[k] -= rt * logStandConc_k;
+    }
+
+    // compute Delta G^0 for all reversible reactions
+    m_rxnstoich->getRevReactionDelta(m_ii, &m_grt[0], &m_rkc[0]);
+ 
+    //doublereal logStandConc = m_kdata->m_logStandConc;
+    doublereal rrt = 1.0/(GasConstant * thermo().temperature());
+    for (i = 0; i < m_nrev; i++) {
+      irxn = m_revindex[i];
+      m_rkc[irxn] = exp(m_rkc[irxn]*rrt);
+    }
+
+    for(i = 0; i != m_nirrev; ++i) {
+      m_rkc[ m_irrev[i] ] = 0.0;
+    }
+  }
+
+  /**
+   * Get the equilibrium constants of all reactions, whether
+   * reversible or not.
+   */
+  void AqueousKinetics::getEquilibriumConstants(doublereal* kc) {
+    int i;
+    _update_rates_T();
+    vector_fp& rkc = m_kdata->m_rkcn;
+  
+    thermo().getStandardChemPotentials(&m_grt[0]);
+    fill(rkc.begin(), rkc.end(), 0.0);
+    doublereal rt = GasConstant * m_kdata->m_temp;
+    int nsp = thermo().nSpecies();
+    for (int k = 0; k < nsp; k++) {
+      doublereal logStandConc_k = thermo().logStandardConc(k);
+      m_grt[k] -= rt * logStandConc_k;
+    }
+    
+    // compute Delta G^0 for all reactions
+    m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], &rkc[0]);
+ 
+    doublereal rrt = 1.0/(GasConstant * thermo().temperature());
+    for (i = 0; i < m_ii; i++) {
+      kc[i] = exp(-rkc[i]*rrt);
+    }
+
+    // force an update of T-dependent properties, so that m_rkcn will
+    // be updated before it is used next.
+    m_kdata->m_temp = 0.0;
+  }
+
+  /**
+   *
+   * getDeltaGibbs():
+   *
+   * Return the vector of values for the reaction gibbs free energy
+   * change
+   * These values depend upon the concentration
+   * of the ideal gas.
+   *
+   *  units = J kmol-1
+   */
+  void AqueousKinetics::getDeltaGibbs(doublereal* deltaG) {
+    /*
+     * Get the chemical potentials of the species in the 
+     * ideal gas solution.
+     */
+    thermo().getChemPotentials(&m_grt[0]);
+    /*
+     * Use the stoichiometric manager to find deltaG for each
+     * reaction.
+     */
+    m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaG);
+  }
+    
+  /**
+   *
+   * getDeltaEnthalpy():
+   * 
+   * Return the vector of values for the reactions change in
+   * enthalpy.
+   * These values depend upon the concentration
+   * of the solution.
+   *
+   *  units = J kmol-1
+   */
+  void AqueousKinetics::getDeltaEnthalpy(doublereal* deltaH) {
+    /*
+     * Get the partial molar enthalpy of all species in the 
+     * ideal gas.
+     */
+    thermo().getPartialMolarEnthalpies(&m_grt[0]);
+    /*
+     * Use the stoichiometric manager to find deltaG for each
+     * reaction.
+     */
+    m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaH);
+  }
+
+  /*
+   *
+   * getDeltaEntropy():
+   *
+   * Return the vector of values for the reactions change in
+   * entropy.
+   * These values depend upon the concentration
+   * of the solution.
+   *
+   *  units = J kmol-1 Kelvin-1
+   */
+  void AqueousKinetics::getDeltaEntropy( doublereal* deltaS) {
+    /*
+     * Get the partial molar entropy of all species in the
+     * solid solution.
+     */
+    thermo().getPartialMolarEntropies(&m_grt[0]);
+    /*
+     * Use the stoichiometric manager to find deltaS for each
+     * reaction.
+     */
+    m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaS);
+  }
+
+  /**
+   *
+   * getDeltaSSGibbs():
+   *
+   * Return the vector of values for the reaction 
+   * standard state gibbs free energy change.
+   * These values don't depend upon the concentration
+   * of the solution.
+   *
+   *  units = J kmol-1
+   */
+  void AqueousKinetics::getDeltaSSGibbs(doublereal* deltaG) {
+    /*
+     *  Get the standard state chemical potentials of the species.
+     *  This is the array of chemical potentials at unit activity 
+     *  We define these here as the chemical potentials of the pure
+     *  species at the temperature and pressure of the solution.
+     */
+    thermo().getStandardChemPotentials(&m_grt[0]);
+    /*
+     * Use the stoichiometric manager to find deltaG for each
+     * reaction.
+     */
+    m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaG);
+  }
+
+  /**
+   *
+   * getDeltaSSEnthalpy():
+   *
+   * Return the vector of values for the change in the
+   * standard state enthalpies of reaction.
+   * These values don't depend upon the concentration
+   * of the solution.
+   *
+   *  units = J kmol-1
+   */
+  void AqueousKinetics::getDeltaSSEnthalpy(doublereal* deltaH) {
+    /*
+     *  Get the standard state enthalpies of the species.
+     *  This is the array of chemical potentials at unit activity 
+     *  We define these here as the enthalpies of the pure
+     *  species at the temperature and pressure of the solution.
+     */
+    thermo().getEnthalpy_RT(&m_grt[0]);
+    doublereal RT = thermo().temperature() * GasConstant;
+    for (int k = 0; k < m_kk; k++) {
+      m_grt[k] *= RT;
+    }
+    /*
+     * Use the stoichiometric manager to find deltaG for each
+     * reaction.
+     */
+    m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaH);
+  }
+
+  /*
+   *
+   * getDeltaSSEntropy():
+   *
+   * Return the vector of values for the change in the
+   * standard state entropies for each reaction.
+   * These values don't depend upon the concentration
+   * of the solution.
+   *
+   *  units = J kmol-1 Kelvin-1
+   */
+  void AqueousKinetics::getDeltaSSEntropy(doublereal* deltaS) {
+    /*
+     *  Get the standard state entropy of the species.
+     *  We define these here as the entropies of the pure
+     *  species at the temperature and pressure of the solution.
+     */
+    thermo().getEntropy_R(&m_grt[0]);
+    doublereal R = GasConstant;
+    for (int k = 0; k < m_kk; k++) {
+      m_grt[k] *= R;
+    }
+    /*
+     * Use the stoichiometric manager to find deltaS for each
+     * reaction.
+     */
+    m_rxnstoich->getReactionDelta(m_ii, &m_grt[0], deltaS);
+  }
+
+
+
+  void AqueousKinetics::updateROP() {
+
+    _update_rates_T();
+    _update_rates_C();
+
+    if (m_kdata->m_ROP_ok) return;
+
+    const vector_fp& rf = m_kdata->m_rfn;
+    const vector_fp& m_rkc = m_kdata->m_rkcn;
+    array_fp& ropf = m_kdata->m_ropf;
+    array_fp& ropr = m_kdata->m_ropr;
+    array_fp& ropnet = m_kdata->m_ropnet;
+
+    // copy rate coefficients into ropf
+    copy(rf.begin(), rf.end(), ropf.begin());
+
+    
+    // multiply by perturbation factor
+    multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
+           
+    // copy the forward rates to the reverse rates                
+    copy(ropf.begin(), ropf.end(), ropr.begin());
+        
+    // for reverse rates computed from thermochemistry, multiply
+    // the forward rates copied into m_ropr by the reciprocals of
+    // the equilibrium constants
+    multiply_each(ropr.begin(), ropr.end(), m_rkc.begin());
+
+    // multiply ropf by concentration products
+    m_rxnstoich->multiplyReactants(&m_conc[0], &ropf[0]); 
+    //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); 
+
+    // for reversible reactions, multiply ropr by concentration
+    // products
+    m_rxnstoich->multiplyRevProducts(&m_conc[0], &ropr[0]); 
+    //m_revProductStoich.multiply(m_conc.begin(), ropr.begin());
+
+    for (int j = 0; j != m_ii; ++j) {
+      ropnet[j] = ropf[j] - ropr[j];
+    }
+
+    m_kdata->m_ROP_ok = true;
+  }
+
+  /**
+   *
+   * getFwdRateConstants():
+   *
+   * Update the rate of progress for the reactions.
+   * This key routine makes sure that the rate of progress vectors
+   * located in the solid kinetics data class are up to date.
+   */
+  void AqueousKinetics::
+  getFwdRateConstants(doublereal *kfwd) {
+    _update_rates_T();
+    _update_rates_C();
+
+    // copy rate coefficients into ropf
+    const vector_fp& rf = m_kdata->m_rfn;
+    array_fp& ropf = m_kdata->m_ropf;
+    copy(rf.begin(), rf.end(), ropf.begin());
+
+   
+
+    // multiply by perturbation factor
+    multiply_each(ropf.begin(), ropf.end(), m_perturb.begin());
+       
+    for (int i = 0; i < m_ii; i++) {
+      kfwd[i] = ropf[i];
+    }
+  }
+
+  /**
+   *
+   * getRevRateConstants():
+   *
+   * Return a vector of the reverse reaction rate constants
+   *
+   * Length is the number of reactions. units depends
+   * on many issues. Note, this routine will return rate constants
+   * for irreversible reactions if the default for
+   * doIrreversible is overridden.
+   */
+  void AqueousKinetics::
+  getRevRateConstants(doublereal *krev, bool doIrreversible) {
+    /*
+     * go get the forward rate constants. -> note, we don't
+     * really care about speed or redundancy in these
+     * informational routines.
+     */
+    getFwdRateConstants(krev);
+
+    if (doIrreversible) {
+      doublereal *tmpKc = &m_kdata->m_ropnet[0];
+      getEquilibriumConstants(tmpKc);
+      for (int i = 0; i < m_ii; i++) {
+	krev[i] /=  tmpKc[i];
+      }
+    } else {
+      /*
+       * m_rkc[] is zero for irreversibly reactions
+       */
+      const vector_fp& m_rkc = m_kdata->m_rkcn;
+      for (int i = 0; i < m_ii; i++) {
+	krev[i] *= m_rkc[i];
+      }
+    }
+  }
+
+  void AqueousKinetics::addReaction(const ReactionData& r) {
+
+    if (r.reactionType == ELEMENTARY_RXN)      addElementaryReaction(r);
+ 
+    // operations common to all reaction types
+    installReagents( r );
+    installGroups(reactionNumber(), r.rgroups, r.pgroups);
+    incrementRxnCount();
+    m_rxneqn.push_back(r.equation);
+  }
+
+
+
+
+  void AqueousKinetics::addElementaryReaction(const ReactionData& r) {
+    int iloc;
+
+    // install rate coeff calculator
+    iloc = m_rates.install( reactionNumber(),
+			    r.rateCoeffType, r.rateCoeffParameters.size(), 
+			    DATA_PTR(r.rateCoeffParameters) );
+
+    // add constant term to rate coeff value vector
+    m_kdata->m_rfn.push_back(r.rateCoeffParameters[0]);                
+
+    // forward rxn order equals number of reactants
+    m_fwdOrder.push_back(r.reactants.size());
+    registerReaction( reactionNumber(), ELEMENTARY_RXN, iloc);
+  }
+
+
+  
+        
+  void AqueousKinetics::installReagents(const ReactionData& r) {
+            
+    m_kdata->m_ropf.push_back(0.0);     // extend by one for new rxn
+    m_kdata->m_ropr.push_back(0.0);
+    m_kdata->m_ropnet.push_back(0.0);
+    int n, ns, m;
+    doublereal nsFlt;
+    doublereal reactantGlobalOrder = 0.0;
+    doublereal productGlobalOrder  = 0.0;
+    int rnum = reactionNumber();
+
+    vector_int rk;
+    int nr = r.reactants.size();
+    for (n = 0; n < nr; n++) {
+      nsFlt = r.rstoich[n];
+      reactantGlobalOrder += nsFlt;
+      ns = (int) nsFlt;
+      if ((doublereal) ns != nsFlt) {
+	if (ns < 1) {
+	  ns = 1;
+	}
+      }
+      if (r.rstoich[n] != 0.0) 
+	m_rrxn[r.reactants[n]][rnum] += r.rstoich[n];
+      for (m = 0; m < ns; m++) {
+	rk.push_back(r.reactants[n]);
+      }
+    }
+    m_reactants.push_back(rk);
+
+    vector_int pk;
+    int np = r.products.size();
+    for (n = 0; n < np; n++) {
+      nsFlt = r.pstoich[n];
+      productGlobalOrder += nsFlt;
+      ns = (int) nsFlt;
+      if ((double) ns != nsFlt) {
+	if (ns < 1) {
+	  ns = 1;
+	}
+      }
+      if (r.pstoich[n] != 0.0) 
+	m_prxn[r.products[n]][rnum] += r.pstoich[n];
+      for (m = 0; m < ns; m++) {
+	pk.push_back(r.products[n]);
+      }
+    }
+    m_products.push_back(pk);
+
+    m_kdata->m_rkcn.push_back(0.0);
+
+    m_rxnstoich->add(reactionNumber(), r);
+
+    if (r.reversible) {
+      m_dn.push_back(productGlobalOrder - reactantGlobalOrder);
+      m_revindex.push_back(reactionNumber());
+      m_nrev++;
+    }
+    else {
+      m_dn.push_back(productGlobalOrder - reactantGlobalOrder);
+      m_irrev.push_back( reactionNumber() );
+      m_nirrev++;
+    }        
+  }
+
+
+  void AqueousKinetics::installGroups(int irxn, 
+				      const vector<grouplist_t>& r, 
+				      const vector<grouplist_t>& p) {
+    if (!r.empty()) {
+      writelog("installing groups for reaction "+int2str(reactionNumber()));
+      m_rgroups[reactionNumber()] = r;
+      m_pgroups[reactionNumber()] = p;
+    }
+  }
+
+
+  void AqueousKinetics::init() { 
+    m_kk = thermo().nSpecies();
+    m_rrxn.resize(m_kk);
+    m_prxn.resize(m_kk);
+    m_conc.resize(m_kk);
+    m_grt.resize(m_kk);
+    m_kdata->m_logp_ref = log(thermo().refPressure()) - log(GasConstant);
+  }
+
+  void AqueousKinetics::finalize() {
+    if (!m_finalized) {
+      m_finalized = true;
+    }
+  }
+
+  bool AqueousKinetics::ready() const {
+    return (m_finalized);
+  }
+
+}

Cantera/src/kinetics/AqueousKinetics.h

@@ -0,0 +1,405 @@
+/**
+ * @file AqueousKinetics.h
+ *
+ * @ingroup chemkinetics
+ *
+ * $Author$
+ * $Revision$
+ * $Date$
+ */
+
+// Copyright 2001  California Institute of Technology
+
+
+#ifndef CT_AQUEOUSKINETICS_H
+#define CT_AQUEOUSKINETICS_H
+
+#include <fstream>
+#include <map>
+
+#include "mix_defs.h"
+#include "Kinetics.h"
+
+#include "utilities.h"
+
+#include "ReactionStoichMgr.h"
+#include "ThirdBodyMgr.h"
+#include "FalloffMgr.h"
+#include "RateCoeffMgr.h"
+
+#include <cmath>
+#include <cstdlib>
+
+void get_wdot(const doublereal* rop, doublereal* wdot);
+
+namespace Cantera {
+
+  // forward references
+
+  
+  class ReactionData;
+  class AqueousKineticsData;
+  class Thermo;
+
+  /**
+   * Holds mechanism-specific data.
+   */
+  class AqueousKineticsData {
+  public:
+    AqueousKineticsData() :
+      m_logp_ref(0.0),
+      m_logc_ref(0.0),
+      m_ROP_ok(false),
+      m_temp(0.0)
+    {}
+    virtual ~AqueousKineticsData(){}
+
+    doublereal m_logp_ref, m_logc_ref;
+    array_fp m_ropf;
+    array_fp m_ropr, m_ropnet;
+    array_fp m_rfn_low, m_rfn_high;
+    bool m_ROP_ok;
+
+    doublereal m_temp;
+    array_fp  m_rfn;
+      
+    array_fp m_rkcn;
+  };
+
+
+  /**
+   * Kinetics manager for elementary aqueous-phase chemistry. This
+   * kinetics manager implements standard mass-action reaction rate
+   * expressions for liquids
+   *
+   *
+   *   Concentration 
+   *
+   * @ingroup kinetics
+   */
+  class AqueousKinetics : public Kinetics {
+
+  public:
+
+    /**
+     * @name Constructors and General Information
+     */
+    //@{
+    /// Constructor.
+    AqueousKinetics(thermo_t* thermo = 0);
+
+    /// Destructor.
+    virtual ~AqueousKinetics();
+
+    virtual int ID() const { return cAqueousKinetics; }
+    virtual int type() const { return cAqueousKinetics; }
+
+    virtual doublereal reactantStoichCoeff(int k, int i) const {
+      return m_rrxn[k][i];
+    }
+
+    virtual doublereal productStoichCoeff(int k, int i) const {
+      return m_prxn[k][i];
+    }
+
+    //@}
+    /**
+     * @name Reaction Rates Of Progress
+     */
+    //@{
+    /**
+     * Forward rates of progress.
+     * Return the forward rates of progress in array fwdROP, which
+     * must be dimensioned at least as large as the total number
+     * of reactions.
+     */
+    virtual void getFwdRatesOfProgress(doublereal* fwdROP) {
+      updateROP();
+      std::copy(m_kdata->m_ropf.begin(), m_kdata->m_ropf.end(), fwdROP);
+    }
+
+    /**
+     * Reverse rates of progress.
+     * Return the reverse rates of progress in array revROP, which
+     * must be dimensioned at least as large as the total number
+     * of reactions.
+     */
+    virtual void getRevRatesOfProgress(doublereal* revROP) {
+      updateROP();
+      std::copy(m_kdata->m_ropr.begin(), m_kdata->m_ropr.end(), revROP);
+    }
+
+    /**
+     * Net rates of progress.  Return the net (forward - reverse)
+     * rates of progress in array netROP, which must be
+     * dimensioned at least as large as the total number of
+     * reactions.
+     */
+    virtual void getNetRatesOfProgress(doublereal* netROP) {
+      updateROP();
+      std::copy(m_kdata->m_ropnet.begin(), m_kdata->m_ropnet.end(), netROP);
+    }
+
+
+    /**
+     * Equilibrium constants. Return the equilibrium constants of
+     * the reactions in concentration units in array kc, which
+     * must be dimensioned at least as large as the total number
+     * of reactions.
+     */
+    virtual void getEquilibriumConstants(doublereal* kc);
+
+    /**
+     * Return the array of values for the reaction gibbs free energy
+     * change.
+     * These values depend on the species concentrations.
+     *
+     *  units = J kmol-1
+     */
+    virtual void getDeltaGibbs( doublereal* deltaG);
+
+    /**
+     * Return the array of values for the reaction enthalpy change.
+     * These values depend upon the species concentrations.
+     *
+     *  units = J kmol-1
+     */
+    virtual void getDeltaEnthalpy( doublereal* deltaH);
+
+    /**
+     * Return the array of values for the reactions change in
+     * entropy.
+     * These values depend upon the concentration
+     * of the solution.
+     *
+     *  units = J kmol-1 Kelvin-1
+     */
+    virtual void getDeltaEntropy(doublereal* deltaS);
+
+    /**
+     * Return the array of values for the reaction
+     * standard state Gibbs free energy change.
+     * These values do not depend on the species 
+     * concentrations.
+     *
+     *  units = J kmol-1
+     */
+    virtual void getDeltaSSGibbs(doublereal* deltaG);
+
+    /**
+     * Return the array of values for the change in the
+     * standard state enthalpies of reaction.
+     * These values do not depend upon the concentration
+     * of the solution.
+     *
+     *  units = J kmol-1
+     */
+    virtual void getDeltaSSEnthalpy(doublereal* deltaH);
+
+    /**
+     * Return the array of values for the change in the
+     * standard state entropies for each reaction.
+     * These values do not depend upon the concentration
+     * of the solution.
+     *
+     *  units = J kmol-1 Kelvin-1
+     */
+    virtual void getDeltaSSEntropy(doublereal* deltaS);
+
+    //@}
+    /**
+     * @name Species Production Rates
+     */
+    //@{
+
+    //! Return the species net production rates
+    /*!
+     * Species net production rates [kmol/m^3/s]. Return the species
+     * net production rates (creation - destruction) in array
+     * wdot, which must be dimensioned at least as large as the
+     * total number of species.
+     *
+     *  @param net  Array of species production rates.
+     *             units kmol m-3 s-1
+     */
+    virtual void getNetProductionRates(doublereal* net) {
+      updateROP();
+      //#ifdef HWMECH
+      //get_wdot(&m_kdata->m_ropnet[0], net);
+      //#else
+      m_rxnstoich->getNetProductionRates(m_kk, 
+					 &m_kdata->m_ropnet[0], net);
+      //#endif
+    }
+
+    /**
+     * Species creation rates [kmol/m^3]. Return the species
+     * creation rates in array cdot, which must be
+     * dimensioned at least as large as the total number of
+     * species.
+     *
+     */
+    virtual void getCreationRates(doublereal* cdot) {
+      updateROP();
+      m_rxnstoich->getCreationRates(m_kk, &m_kdata->m_ropf[0],
+				    &m_kdata->m_ropr[0], cdot);
+    }
+
+    /**
+     * Species destruction rates [kmol/m^3]. Return the species
+     * destruction rates in array ddot, which must be
+     * dimensioned at least as large as the total number of
+     * species.
+     *
+     */
+    virtual void getDestructionRates(doublereal* ddot) {
+      updateROP();
+      m_rxnstoich->getDestructionRates(m_kk, &m_kdata->m_ropf[0],
+				       &m_kdata->m_ropr[0], ddot);
+           
+    }
+
+    //@}
+    /**
+     * @name Reaction Mechanism Informational Query Routines
+     */
+    //@{
+
+    /**
+     * Flag specifying the type of reaction. The legal values and
+     * their meaning are specific to the particular kinetics
+     * manager.
+     */
+    virtual int reactionType(int i) const {
+      return m_index[i].first;
+    }
+
+    virtual std::string reactionString(int i) const {
+      return m_rxneqn[i];
+    }
+
+    /**
+     * True if reaction i has been declared to be reversible. If
+     * isReversible(i) is false, then the reverse rate of progress
+     * for reaction i is always zero.
+     */
+    virtual bool isReversible(int i) {
+      if (std::find(m_revindex.begin(), m_revindex.end(), i)
+	  < m_revindex.end()) return true;
+      else return false;
+    }
+
+    /**
+     * Return the forward rate constants
+     *
+     * length is the number of reactions. units depends
+     * on many issues.
+     */
+    virtual void getFwdRateConstants(doublereal *kfwd);
+
+    /**
+     * Return the reverse rate constants.
+     *
+     * length is the number of reactions. units depends
+     * on many issues. Note, this routine will return rate constants
+     * for irreversible reactions if the default for
+     * doIrreversible is overridden.
+     */
+    virtual void getRevRateConstants(doublereal *krev,
+				     bool doIrreversible = false);
+
+    //@}
+    /**
+     * @name Reaction Mechanism Setup Routines
+     */
+    //@{
+
+    virtual void init();
+
+    ///  Add a reaction to the mechanism.
+    virtual void addReaction(const ReactionData& r);
+
+    virtual void finalize();
+    virtual bool ready() const;
+
+    virtual void update_T();
+    virtual void update_C();
+
+    void updateROP();
+
+
+    const std::vector<grouplist_t>& reactantGroups(int i)
+    { return m_rgroups[i]; }
+    const std::vector<grouplist_t>& productGroups(int i)
+    { return m_pgroups[i]; }
+
+
+    void _update_rates_T();
+    void _update_rates_C();
+
+    //@}
+
+  protected:
+
+    int                                 m_kk, m_nfall;
+  
+    Rate1<Arrhenius>                    m_rates;
+
+    mutable std::map<int, std::pair<int, int> >   m_index;
+
+    std::vector<int> m_irrev;
+
+    ReactionStoichMgr*                   m_rxnstoich;
+
+    std::vector<int>                         m_fwdOrder;
+
+    int m_nirrev;
+    int m_nrev;
+
+    std::map<int, std::vector<grouplist_t> >      m_rgroups;
+    std::map<int, std::vector<grouplist_t> >      m_pgroups;
+
+    std::vector<int>                         m_rxntype;
+
+    mutable std::vector<std::map<int, doublereal> >     m_rrxn;
+    mutable std::vector<std::map<int, doublereal> >     m_prxn;
+
+    /**
+     * Difference between the input global reactants order
+     * and the input global products order. Changed to a double
+     * to account for the fact that we can have real-valued
+     * stoichiometries.
+     */
+    array_fp  m_dn;
+    array_int m_revindex;
+
+    std::vector<std::string> m_rxneqn;
+
+    AqueousKineticsData* m_kdata;
+
+    array_fp m_conc;
+    array_fp m_grt;
+
+
+  private:
+
+    int reactionNumber(){ return m_ii;}
+    std::vector<std::map<int, doublereal> > m_stoich;
+
+    void addElementaryReaction(const ReactionData& r);
+ 
+
+    void installReagents(const ReactionData& r);
+
+    void installGroups(int irxn, const std::vector<grouplist_t>& r,
+		       const std::vector<grouplist_t>& p);
+    void updateKc();
+
+    void registerReaction(int rxnNumber, int type, int loc) {
+      m_index[rxnNumber] = std::pair<int, int>(type, loc);
+    }
+    bool m_finalized;
+  };
+}
+
+#endif

Cantera/src/kinetics/InterfaceKinetics.cpp

@@ -191,78 +191,78 @@ namespace Cantera {
     }
 
 
-    /**
+  /**
-     * Update properties that depend on concentrations. This method
+   * Update properties that depend on concentrations. This method
-     * fills out the array of generalized concentrations by calling
+   * fills out the array of generalized concentrations by calling
-     * method getActivityConcentrations for each phase, which classes
+   * method getActivityConcentrations for each phase, which classes
-     * representing phases should overload to return the appropriate
+   * representing phases should overload to return the appropriate
-     * quantities.
+   * quantities.
-     */ 
+   */ 
-    void InterfaceKinetics::
+  void InterfaceKinetics::
-    _update_rates_C() {
+  _update_rates_C() {
-        int n;
+    int n;
 
-        int np = nPhases();
+    int np = nPhases();
-        for (n = 0; n < np; n++) {
+    for (n = 0; n < np; n++) {
-	  /*
+      /*
-	   * We call the getActivityConcentrations function of each
+       * We call the getActivityConcentrations function of each
-	   * ThermoPhase class that makes up this kinetics object to 
+       * ThermoPhase class that makes up this kinetics object to 
-	   * obtain the generalized concentrations for species within that 
+       * obtain the generalized concentrations for species within that 
-	   * class. This is collected in the vector m_conc. m_start[]
+       * class. This is collected in the vector m_conc. m_start[]
-	   * are integer indecises for that vector denoting the start of the
+       * are integer indecises for that vector denoting the start of the
-	   * species for each phase.
+       * species for each phase.
-	   */
+       */
-            thermo(n).getActivityConcentrations(DATA_PTR(m_conc) + m_start[n]);
+      thermo(n).getActivityConcentrations(DATA_PTR(m_conc) + m_start[n]);
-        }
-        m_kdata->m_ROP_ok = false;
     }
+    m_kdata->m_ROP_ok = false;
+  }
 
 
-    /**
+  /**
-     * Update the equilibrium constants in molar units for all
+   * Update the equilibrium constants in molar units for all
-     * reversible reactions. Irreversible reactions have their 
+   * reversible reactions. Irreversible reactions have their 
-     * equilibrium constant set to zero.
+   * equilibrium constant set to zero.
-     */
+   */
-    void InterfaceKinetics::updateKc() {
+  void InterfaceKinetics::updateKc() {
-        int i, irxn;
+    int i, irxn;
 
-        vector_fp& m_rkc = m_kdata->m_rkcn;
+    vector_fp& m_rkc = m_kdata->m_rkcn;
-        fill(m_rkc.begin(), m_rkc.end(), 0.0);
+    fill(m_rkc.begin(), m_rkc.end(), 0.0);
 
-        //static vector_fp mu(nTotalSpecies());
+    //static vector_fp mu(nTotalSpecies());
-        if (m_nrev > 0) {
+    if (m_nrev > 0) {
 
-            int n, nsp, k, ik=0;
+      int n, nsp, k, ik = 0;
-            doublereal rt = GasConstant*thermo(0).temperature();
+      doublereal rt = GasConstant*thermo(0).temperature();
-            doublereal rrt = 1.0/rt;
+      doublereal rrt = 1.0 / rt;
-            int np = nPhases();
+      int np = nPhases();
-            for (n = 0; n < np; n++) {
+      for (n = 0; n < np; n++) {
-                thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]);
+	thermo(n).getStandardChemPotentials(DATA_PTR(m_mu0) + m_start[n]);
-                nsp = thermo(n).nSpecies();
+	nsp = thermo(n).nSpecies();
-                for (k = 0; k < nsp; k++) {
+	for (k = 0; k < nsp; k++) {
-                    m_mu0[ik] -= rt*thermo(n).logStandardConc(k);
+	  m_mu0[ik] -= rt * thermo(n).logStandardConc(k);
-                    m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
+	  m_mu0[ik] += Faraday * m_phi[n] * thermo(n).charge(k);
-                    ik++;
+	  ik++;
-                }
+	}
-            }
+      }
 
-            // compute Delta mu^0 for all reversible reactions
+      // compute Delta mu^0 for all reversible reactions
-            m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0), 
+      m_rxnstoich.getRevReactionDelta(m_ii, DATA_PTR(m_mu0), 
-                DATA_PTR(m_rkc));
+				      DATA_PTR(m_rkc));
 
-            for (i = 0; i < m_nrev; i++) {
+      for (i = 0; i < m_nrev; i++) {
-                irxn = m_revindex[i];
+	irxn = m_revindex[i];
-                if (irxn < 0 || irxn >= nReactions()) {
+	if (irxn < 0 || irxn >= nReactions()) {
-                    throw CanteraError("InterfaceKinetics",
+	  throw CanteraError("InterfaceKinetics",
-                        "illegal value: irxn = "+int2str(irxn));
+			     "illegal value: irxn = "+int2str(irxn));
-                }
+	}
-                m_rkc[irxn] = exp(m_rkc[irxn]*rrt);
+	m_rkc[irxn] = exp(m_rkc[irxn]*rrt);
-            }
+      }
-            for (i = 0; i != m_nirrev; ++i) {
+      for (i = 0; i != m_nirrev; ++i) {
-                m_rkc[ m_irrev[i] ] = 0.0;
+	m_rkc[ m_irrev[i] ] = 0.0;
-            }
+      }
-        }
     }
+  }
 
 
     void InterfaceKinetics::checkPartialEquil() {

Cantera/src/kinetics/KineticsFactory.cpp

@@ -22,6 +22,7 @@
 #include "InterfaceKinetics.h"
 #include "EdgeKinetics.h"
 #include "importKinetics.h"
+#include "AqueousKinetics.h"
 
 using namespace std;
 
@@ -33,8 +34,8 @@ namespace Cantera {
       #endif
 
     static int ntypes = 5;
-    static string _types[] = {"none", "GasKinetics", "GRI30", "Interface", "Edge"};
+    static string _types[] = {"none", "GasKinetics", "GRI30", "Interface", "Edge", "AqueousKinetics"};
-    static int _itypes[]   = {0, cGasKinetics, cGRI30, cInterfaceKinetics, cEdgeKinetics};
+    static int _itypes[]   = {0, cGasKinetics, cGRI30, cInterfaceKinetics, cEdgeKinetics, cAqueousKinetics};
 
     /**
      * Return a new kinetics manager that implements a reaction
@@ -106,6 +107,10 @@ namespace Cantera {
             k = new EdgeKinetics;
             break;
       
+         case cAqueousKinetics:
+            k = new AqueousKinetics;
+            break;
+
         default:
         throw UnknownKineticsModel("KineticsFactory::newKinetics", 
                        kintype);

Cantera/src/kinetics/KineticsFactory.h

@@ -31,7 +31,7 @@ namespace Cantera {
         CanteraError(proc, "Specified Kinetics model "   
              + kineticsModel + 
              " does not match any known type.") {}
-    virtual ~UnknownKineticsModel() {}
+    virtual ~UnknownKineticsModel() throw() {}
     };
 
 

Cantera/src/kinetics/Makefile.in

@@ -36,13 +36,14 @@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG)
 
 ifeq ($(do_kinetics),1)
 KINETICS_OBJ=importKinetics.o GRI_30_Kinetics.o KineticsFactory.o \
-             GasKinetics.o \
+             GasKinetics.o AqueousKinetics.o \
              FalloffFactory.o ReactionStoichMgr.o Kinetics.o solveSP.o
 KINETICS_H = importKinetics.h GRI_30_Kinetics.h KineticsFactory.h \
              Kinetics.h GasKinetics.h \
              FalloffFactory.h  ReactionStoichMgr.h reaction_defs.h \
              FalloffMgr.h ThirdBodyMgr.h RateCoeffMgr.h ReactionData.h \
-             RxnRates.h Enhanced3BConc.h StoichManager.h solveSP.h
+             RxnRates.h Enhanced3BConc.h StoichManager.h solveSP.h \
+             AqueousKinetics.h
 KINETICS   = $(KINETICS_OBJ) $(KINETICS_H)
 endif