diff --git a/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc
new file mode 100644
index 00000000..d1e791c1
--- /dev/null
+++ b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc
@@ -0,0 +1,85 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#include "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h"
+
+#include <boost/property_tree/json_parser.hpp>
+#include "Database2JSON.h"
+namespace pt = boost::property_tree;
+
+#include "SAMRAI/tbox/InputManager.h"
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/math/HierarchyCellDataOpsReal.h"
+
+using namespace SAMRAI;
+
+#include <cassert>
+
+//=======================================================================
+
+CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::
+    CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB(
+        const int norderp_A, boost::property_tree::ptree calphad_pt,
+        std::shared_ptr<tbox::Database> newton_db,
+        const Thermo4PFM::ConcInterpolationType conc_interp_func_type,
+        MolarVolumeStrategy* mvstrategy, const int conc_l_id,
+        const int conc_a_id, const int conc_b_id, const int conc_id)
+    : CALPHADFreeEnergyBinaryMultiOrderThreePhases<
+          Thermo4PFM::CALPHADFreeEnergyFunctionsBinaryThreePhase>(
+          calphad_pt, newton_db, conc_interp_func_type, norderp_A, mvstrategy,
+          conc_l_id, conc_a_id, conc_b_id)
+{
+   tbox::plog << "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB..."
+              << std::endl;
+   setup(calphad_pt, newton_db);
+
+   d_multiorder_driving_force.reset(
+       new MultiOrderBinaryThreePhasesDrivingForceStochioAB(this, conc_id,
+                                                            norderp_A));
+}
+
+//=======================================================================
+
+void CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::setup(
+    pt::ptree calphad_pt, std::shared_ptr<tbox::Database> newton_db)
+{
+   pt::ptree newton_pt;
+   if (newton_db) copyDatabase(newton_db, newton_pt);
+   // set looser tol in solver
+   newton_pt.put("tol", 1.e-6);
+   d_ceq_fenergy.reset(new Thermo4PFM::CALPHADFreeEnergyFunctionsBinary(
+       calphad_pt, newton_pt, d_energy_interp_func_type,
+       d_conc_interp_func_type));
+}
+
+
+bool CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::computeCeqT(
+    const double temperature, const Thermo4PFM::PhaseIndex pi0,
+    const Thermo4PFM::PhaseIndex pi1, double* ceq)
+{
+   std::cout << "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::"
+                "computeCeqT..."
+             << std::endl;
+   return d_ceq_fenergy->computeCeqT(temperature, &ceq[0], 50, true);
+}
+
+void CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB::addDrivingForce(
+    const double time, hier::Patch& patch, const int temperature_id,
+    const int phase_id, const int conc_id, const int f_l_id, const int f_a_id,
+    const int f_b_id, const int rhs_id)
+{
+   (void)time;
+
+   d_multiorder_driving_force->addDrivingForce(patch, temperature_id, phase_id,
+                                               d_conc_l_id, d_conc_a_id,
+                                               d_conc_b_id, f_l_id, f_a_id,
+                                               f_b_id, rhs_id);
+}
diff --git a/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h
new file mode 100644
index 00000000..8b4d7256
--- /dev/null
+++ b/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h
@@ -0,0 +1,61 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#ifndef included_CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB
+#define included_CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB
+
+#include "CALPHADFreeEnergyBinaryMultiOrderThreePhases.h"
+#include "InterpolationType.h"
+#include "MultiOrderBinaryThreePhasesDrivingForceStochioAB.h"
+
+#include "CALPHADFreeEnergyFunctionsBinaryThreePhase.h"
+#include "CALPHADFreeEnergyFunctionsBinary.h"
+
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/tbox/Database.h"
+#include "SAMRAI/hier/Box.h"
+class MolarVolumeStrategy;
+
+#include <string>
+#include <vector>
+#include <boost/property_tree/ptree.hpp>
+
+class CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB
+    : public CALPHADFreeEnergyBinaryMultiOrderThreePhases<
+          Thermo4PFM::CALPHADFreeEnergyFunctionsBinaryThreePhase>
+{
+ public:
+   CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB(
+       const int norderp_A, boost::property_tree::ptree calphad_db,
+       std::shared_ptr<tbox::Database> newton_db,
+       const Thermo4PFM::ConcInterpolationType conc_interp_func_type,
+       MolarVolumeStrategy* mvstrategy, const int conc_l_id,
+       const int conc_a_id, const int conc_b_id, const int conc_id);
+
+   bool computeCeqT(const double temperature, const Thermo4PFM::PhaseIndex pi0,
+                    const Thermo4PFM::PhaseIndex pi1, double* ceq) override;
+
+   void addDrivingForce(const double time, hier::Patch& patch,
+                        const int temperature_id, const int phase_id,
+                        const int conc_id, const int f_l_id, const int f_a_id,
+                        const int f_b_id, const int rhs_id) override;
+
+ private:
+   void setup(boost::property_tree::ptree calphad_pt,
+              std::shared_ptr<tbox::Database> newton_db);
+
+   std::shared_ptr<Thermo4PFM::CALPHADFreeEnergyFunctionsBinary> d_ceq_fenergy;
+
+
+   std::shared_ptr<MultiOrderBinaryThreePhasesDrivingForceStochioAB>
+       d_multiorder_driving_force;
+};
+
+#endif
diff --git a/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc
new file mode 100644
index 00000000..6d21f4e5
--- /dev/null
+++ b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc
@@ -0,0 +1,87 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#include "Database2JSON.h"
+
+#include <boost/property_tree/json_parser.hpp>
+namespace pt = boost::property_tree;
+
+#include "CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h"
+
+CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB ::
+    CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB(
+        const short norderp_A, const int conc_l_id, const int conc_a_id,
+        const int conc_b_id, const QuatModelParameters& model_parameters,
+        std::shared_ptr<tbox::Database> conc_db)
+    : EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases(
+          norderp_A, conc_l_id, conc_a_id, conc_b_id, model_parameters,
+          conc_db),
+      d_model_parameters(model_parameters)
+{
+   tbox::plog << "CALPHADequilPhaseConcMultiOrderThreePhasesStochioAB..."
+              << std::endl;
+   std::shared_ptr<tbox::Database> conc_calphad_db =
+       conc_db->getDatabase("Calphad");
+   std::string calphad_filename = conc_calphad_db->getString("filename");
+
+   std::shared_ptr<tbox::MemoryDatabase> calphad_db;
+   boost::property_tree::ptree calphad_pt;
+
+   if (calphad_filename.compare(calphad_filename.size() - 4, 4, "json") == 0) {
+      boost::property_tree::read_json(calphad_filename, calphad_pt);
+   } else {
+      calphad_db.reset(new tbox::MemoryDatabase("calphad_db"));
+      tbox::InputManager::getManager()->parseInputFile(calphad_filename,
+                                                       calphad_db);
+      copyDatabase(calphad_db, calphad_pt);
+   }
+}
+
+int CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB ::
+    computePhaseConcentrations(const double temp, double* c, double* hphi,
+                               double* x)
+{
+   assert(!std::isnan(x[0]));
+
+   const double epsilon = 1e-1;
+
+   const double cA = d_model_parameters.getStochioA();
+   const double cB = d_model_parameters.getStochioB();
+
+   const double xeq = d_model_parameters.ceq_liquid(temp);
+   // std::cout << "xeq = " << xeq << std::endl;
+   assert(!std::isnan(xeq));
+
+   // solve explicit KKS problem
+   const double a = hphi[0] + epsilon * std::exp(-hphi[0] / epsilon);
+   const double dkks = (c[0] - hphi[0] * xeq - hphi[1] * cA - hphi[2] * cB) / a;
+
+   const double xmin = 0.7857;
+   const double xmax = 0.999;
+   const double d = dkks > 0. ? (xmax - xeq) : xeq - xmin;
+   const double t = std::tanh(dkks / d);
+   const double xkks = xeq + d * t;
+
+#if 0
+   for (short i = 0; i < 3; i++)
+      std::cerr << hphi[i] << ", ";
+   std::cerr << ", c=" << c[0] << std::endl;
+   std::cerr << "x = " << x[0] << std::endl;
+   std::cerr << "xkks = " << xkks << std::endl;
+   std::cerr << "conc[0] - hphi1 * cA - hphi2 * cB = "
+             << c[0] - hphi[1] * cA - hphi[2] * cB << std::endl;
+#endif
+
+   x[0] = xkks;
+   x[1] = cA;
+   x[2] = cB;
+
+   return 0;
+}
diff --git a/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h
new file mode 100644
index 00000000..83a2649a
--- /dev/null
+++ b/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h
@@ -0,0 +1,37 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#ifndef included_CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB
+#define included_CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB
+
+#include "EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.h"
+#include "InterpolationType.h"
+#include "QuatModelParameters.h"
+
+#include "SAMRAI/tbox/InputManager.h"
+
+class CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB
+    : public EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases
+{
+ public:
+   CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB(
+       const short norderp_A, const int conc_l_id, const int conc_a_id,
+       const int conc_b_id, const QuatModelParameters& model_parameters,
+       std::shared_ptr<tbox::Database> conc_db);
+
+   ~CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB() {}
+
+ private:
+   QuatModelParameters d_model_parameters;
+
+   int computePhaseConcentrations(double, double*, double*, double*) override;
+};
+
+#endif
diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt
index 145ec860..cd468efb 100644
--- a/source/CMakeLists.txt
+++ b/source/CMakeLists.txt
@@ -30,6 +30,7 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc
             ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryDrivingForce.cc
             ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryThreePhasesDrivingForce.cc
 	    ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryThreePhasesDrivingForceStochioB.cc
+	    ${CMAKE_SOURCE_DIR}/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc
             ${CMAKE_SOURCE_DIR}/source/CahnHilliardDoubleWell.cc
             ${CMAKE_SOURCE_DIR}/source/SinteringUWangRHSStrategy.cc
             ${CMAKE_SOURCE_DIR}/source/WangSinteringCompositionRHSStrategy.cc
@@ -79,12 +80,14 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc
             ${CMAKE_SOURCE_DIR}/source/QuadraticFreeEnergyMultiOrderBinaryThreePhase.cc
             ${CMAKE_SOURCE_DIR}/source/QuadraticFreeEnergyMultiOrderTernaryThreePhase.cc
             ${CMAKE_SOURCE_DIR}/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc
+            ${CMAKE_SOURCE_DIR}/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.cc
             ${CMAKE_SOURCE_DIR}/source/EquilibriumPhaseConcentrationsBinaryMultiOrder.cc
             ${CMAKE_SOURCE_DIR}/source/EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.cc
             ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsStrategyMultiOrder.cc
             ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsStrategyMultiOrderThreePhases.cc
             ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsThreePhasesStochioB.cc
             ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioB.cc
+            ${CMAKE_SOURCE_DIR}/source/CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.cc
             ${CMAKE_SOURCE_DIR}/source/QuadraticEquilibriumPhaseConcentrationsBinaryMultiOrder.cc
             ${CMAKE_SOURCE_DIR}/source/QuadraticEquilibriumThreePhasesTernaryMultiOrder.cc
             ${CMAKE_SOURCE_DIR}/source/QuadraticEquilibriumPhaseConcentrationsBinary.cc
@@ -93,7 +96,9 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc
             ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyStrategyBinaryFolchPlapp.cc
             ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyStrategyBinaryFolchPlappStochioB.cc
             ${CMAKE_SOURCE_DIR}/source/ParabolicEquilibriumThreePhasesBinaryMultiOrder.cc
+            ${CMAKE_SOURCE_DIR}/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.cc
             ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioB.cc
+            ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.cc
             ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyStrategyTernary.cc
             ${CMAKE_SOURCE_DIR}/source/tools.cc
             ${CMAKE_SOURCE_DIR}/source/FreeEnergyStrategy.cc
diff --git a/source/FreeEnergyStrategyFactory.h b/source/FreeEnergyStrategyFactory.h
index 03e3dd3e..6fb29675 100644
--- a/source/FreeEnergyStrategyFactory.h
+++ b/source/FreeEnergyStrategyFactory.h
@@ -20,6 +20,7 @@
 #include "CALPHADFreeEnergyStrategyBinaryFolchPlapp.h"
 #include "CALPHADFreeEnergyBinaryMultiOrderThreePhases.h"
 #include "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioB.h"
+#include "CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB.h"
 #include "CALPHADFreeEnergyStrategyBinaryFolchPlappStochioB.h"
 #include "ParabolicFreeEnergyBinary.h"
 #include "ParabolicFreeEnergyMultiOrderBinary.h"
@@ -27,6 +28,7 @@
 #include "QuadraticFreeEnergyMultiOrderBinary.h"
 #include "QuadraticFreeEnergyMultiOrderTernaryThreePhase.h"
 #include "ParabolicFreeEnergyMultiOrderBinaryThreePhase.h"
+#include "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h"
 #include "KKSdiluteBinary.h"
 #include "BiasDoubleWellBeckermannFreeEnergyStrategy.h"
 #include "BiasDoubleWellUTRCFreeEnergyStrategy.h"
@@ -41,7 +43,8 @@ class FreeEnergyStrategyFactory
    static std::shared_ptr<FreeEnergyStrategy> create(
        QuatModelParameters& model_parameters, const int ncompositions,
        const int conc_l_scratch_id, const int conc_a_scratch_id,
-       const int conc_b_scratch_id, MolarVolumeStrategy* mvstrategy,
+       const int conc_b_scratch_id, const int conc_scratch_id,
+       MolarVolumeStrategy* mvstrategy,
        MeltingTemperatureStrategy* meltingT_strategy, const double Tref,
        std::shared_ptr<tbox::Database> conc_db)
    {
@@ -51,6 +54,8 @@ class FreeEnergyStrategyFactory
       std::shared_ptr<FreeEnergyStrategy> free_energy_strategy;
 
       if (model_parameters.with_concentration()) {
+         const double cB = model_parameters.getStochioB();
+         const double cA = model_parameters.getStochioA();
 
          if (model_parameters.isConcentrationModelCALPHAD()) {
             std::shared_ptr<tbox::MemoryDatabase> calphad_db;
@@ -82,7 +87,17 @@ class FreeEnergyStrategyFactory
                if (model_parameters.withMultipleOrderP()) {
                   tbox::plog << "MultiOrder..." << std::endl;
                   if (conc_b_scratch_id >= 0) {
-                     if (model_parameters.getStochioB() >= 0.) {
+                     if (cA >= 0. && cB >= 0.) {
+                        tbox::plog << "StochioAB..." << std::endl;
+                        free_energy_strategy.reset(
+                            new CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioAB(
+                                model_parameters.norderpA(), calphad_pt,
+                                newton_db,
+                                model_parameters.conc_interp_func_type(),
+                                mvstrategy, conc_l_scratch_id,
+                                conc_a_scratch_id, conc_b_scratch_id,
+                                conc_scratch_id));
+                     } else if (model_parameters.getStochioB() >= 0.) {
                         tbox::plog << "StochioB..." << std::endl;
                         free_energy_strategy.reset(
                             new CALPHADFreeEnergyBinaryMultiOrderThreePhasesStochioB(
@@ -208,13 +223,22 @@ class FreeEnergyStrategyFactory
                   tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhase."
                                 ".."
                              << std::endl;
-                  free_energy_strategy.reset(
-                      new ParabolicFreeEnergyMultiOrderBinaryThreePhase(
-                          conc_db->getDatabase("Parabolic"),
-                          model_parameters.conc_interp_func_type(),
-                          model_parameters.norderpA(), mvstrategy,
-                          conc_l_scratch_id, conc_a_scratch_id,
-                          conc_b_scratch_id));
+                  if (cA > 0. && cB > 0.) {
+                     free_energy_strategy.reset(
+                         new ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB(
+                             conc_db->getDatabase("Parabolic"),
+                             model_parameters.norderpA(), mvstrategy,
+                             conc_l_scratch_id, conc_a_scratch_id,
+                             conc_b_scratch_id));
+
+                  } else {
+                     free_energy_strategy.reset(
+                         new ParabolicFreeEnergyMultiOrderBinaryThreePhase(
+                             conc_db->getDatabase("Parabolic"),
+                             model_parameters.norderpA(), mvstrategy,
+                             conc_l_scratch_id, conc_a_scratch_id,
+                             conc_b_scratch_id));
+                  }
                } else {
                   tbox::plog << "ParabolicFreeEnergyMultiOrderBinary"
                              << std::endl;
diff --git a/source/MultiOrderBinaryDrivingForce.h b/source/MultiOrderBinaryDrivingForce.h
index 4f5fd0ab..c1960256 100644
--- a/source/MultiOrderBinaryDrivingForce.h
+++ b/source/MultiOrderBinaryDrivingForce.h
@@ -15,14 +15,12 @@
 
 #include "SAMRAI/hier/Patch.h"
 
-using namespace SAMRAI;
-
 class MultiOrderBinaryDrivingForce
 {
  public:
    MultiOrderBinaryDrivingForce(FreeEnergyStrategyBinary* fenergy_strategy);
 
-   void addDrivingForce(hier::Patch& patch, const int temperature_id,
+   void addDrivingForce(SAMRAI::hier::Patch& patch, const int temperature_id,
                         const int phase_id, const int conc_l_id,
                         const int conc_a_id, const int f_l_id, const int f_a_id,
                         const int rhs_id);
diff --git a/source/MultiOrderBinaryThreePhasesDrivingForce.h b/source/MultiOrderBinaryThreePhasesDrivingForce.h
index 294f77b6..aaa02e29 100644
--- a/source/MultiOrderBinaryThreePhasesDrivingForce.h
+++ b/source/MultiOrderBinaryThreePhasesDrivingForce.h
@@ -15,15 +15,13 @@
 
 #include "SAMRAI/hier/Patch.h"
 
-using namespace SAMRAI;
-
 class MultiOrderBinaryThreePhasesDrivingForce
 {
  public:
    MultiOrderBinaryThreePhasesDrivingForce(
        FreeEnergyStrategyBinary* fenergy_strategy, const int norderp_A);
 
-   void addDrivingForce(hier::Patch& patch, const int temperature_id,
+   void addDrivingForce(SAMRAI::hier::Patch& patch, const int temperature_id,
                         const int phase_id, const int conc_l_id,
                         const int conc_a_id, const int conc_b_id,
                         const int f_l_id, const int f_a_id, const int f_b_id,
diff --git a/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc
new file mode 100644
index 00000000..334e0250
--- /dev/null
+++ b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.cc
@@ -0,0 +1,290 @@
+#include "MultiOrderBinaryThreePhasesDrivingForceStochioAB.h"
+
+#include "functions.h"
+#include "SAMRAI/pdat/CellData.h"
+
+
+using namespace SAMRAI;
+
+MultiOrderBinaryThreePhasesDrivingForceStochioAB::
+    MultiOrderBinaryThreePhasesDrivingForceStochioAB(
+        FreeEnergyStrategyBinary* fenergy_strategy, const int conc_id,
+        const int norderp_A)
+    : d_fenergy_strategy(fenergy_strategy),
+      d_conc_id(conc_id),
+      d_norderp_A(norderp_A)
+{
+}
+
+void MultiOrderBinaryThreePhasesDrivingForceStochioAB::addDrivingForce(
+    hier::Patch& patch, const int temperature_id, const int phase_id,
+    const int conc_l_id, const int conc_a_id, const int conc_b_id,
+    const int f_l_id, const int f_a_id, const int f_b_id, const int rhs_id)
+{
+   assert(phase_id >= 0);
+   assert(f_l_id >= 0);
+   assert(f_a_id >= 0);
+   assert(f_b_id >= 0);
+   assert(rhs_id >= 0);
+   assert(conc_l_id >= 0);
+   assert(conc_a_id >= 0);
+   assert(conc_b_id >= 0);
+   assert(temperature_id >= 0);
+   assert(d_conc_id >= 0);
+
+   std::shared_ptr<pdat::CellData<double> > phase(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(phase_id)));
+   assert(phase);
+   assert(phase->getDepth() > 1);
+
+   std::shared_ptr<pdat::CellData<double> > conc(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(d_conc_id)));
+   assert(conc);
+
+   std::shared_ptr<pdat::CellData<double> > temperature(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(temperature_id)));
+   assert(temperature);
+
+   std::shared_ptr<pdat::CellData<double> > fl(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(f_l_id)));
+   assert(fl);
+
+   std::shared_ptr<pdat::CellData<double> > fa(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(f_a_id)));
+   assert(fa);
+
+   std::shared_ptr<pdat::CellData<double> > fb(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(f_b_id)));
+   assert(fb);
+
+   std::shared_ptr<pdat::CellData<double> > cl(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(conc_l_id)));
+   assert(cl);
+
+   std::shared_ptr<pdat::CellData<double> > ca(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(conc_a_id)));
+   assert(ca);
+
+   std::shared_ptr<pdat::CellData<double> > cb(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(conc_b_id)));
+   assert(cb);
+
+   std::shared_ptr<pdat::CellData<double> > rhs(
+       SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+           patch.getPatchData(rhs_id)));
+
+   assert(fl->getGhostCellWidth()[0] == fa->getGhostCellWidth()[0]);
+   assert(cl->getGhostCellWidth()[0] == ca->getGhostCellWidth()[0]);
+   if (conc->getGhostCellWidth()[0] != phase->getGhostCellWidth()[0])
+      std::cerr << "Wrong number of ghosts!!!" << std::endl;
+   assert(phase->getDepth() > 1);
+   assert(phase->getDepth() == rhs->getDepth());
+
+   const hier::Box& pbox(patch.getBox());
+   const int norderp = phase->getDepth();
+
+   double* ptr_temp = temperature->getPointer();
+   double* ptr_f_l = fl->getPointer();
+   double* ptr_f_a = fa->getPointer();
+   double* ptr_f_b = fb->getPointer();
+   double* ptr_c_l = cl->getPointer();
+   double* ptr_c_a = ca->getPointer();
+   double* ptr_c_b = cb->getPointer();
+   double* ptr_c = conc->getPointer();
+
+   const hier::Box& rhs_gbox = rhs->getGhostBox();
+   int imin_rhs = rhs_gbox.lower(0);
+   int jmin_rhs = rhs_gbox.lower(1);
+   int jp_rhs = rhs_gbox.numberCells(0);
+   int kmin_rhs = 0;
+   int kp_rhs = 0;
+#if (NDIM == 3)
+   kmin_rhs = rhs_gbox.lower(2);
+   kp_rhs = jp_rhs * rhs_gbox.numberCells(1);
+#endif
+
+   const hier::Box& temp_gbox = temperature->getGhostBox();
+   int imin_temp = temp_gbox.lower(0);
+   int jmin_temp = temp_gbox.lower(1);
+   int jp_temp = temp_gbox.numberCells(0);
+   int kmin_temp = 0;
+   int kp_temp = 0;
+#if (NDIM == 3)
+   kmin_temp = temp_gbox.lower(2);
+   kp_temp = jp_temp * temp_gbox.numberCells(1);
+#endif
+
+   const hier::Box& pf_gbox = phase->getGhostBox();
+   int imin_pf = pf_gbox.lower(0);
+   int jmin_pf = pf_gbox.lower(1);
+   int jp_pf = pf_gbox.numberCells(0);
+   int kmin_pf = 0;
+   int kp_pf = 0;
+#if (NDIM == 3)
+   kmin_pf = pf_gbox.lower(2);
+   kp_pf = jp_pf * pf_gbox.numberCells(1);
+#endif
+
+   // Assuming f_l, f_a, all have same ghost box
+   const hier::Box& f_i_gbox = fl->getGhostBox();
+   int imin_f_i = f_i_gbox.lower(0);
+   int jmin_f_i = f_i_gbox.lower(1);
+   int jp_f_i = f_i_gbox.numberCells(0);
+   int kmin_f_i = 0;
+   int kp_f_i = 0;
+#if (NDIM == 3)
+   kmin_f_i = f_i_gbox.lower(2);
+   kp_f_i = jp_f_i * f_i_gbox.numberCells(1);
+#endif
+
+   // Assuming c_l, c_a, all have same ghost box
+   const hier::Box& c_i_gbox = cl->getGhostBox();
+   int imin_c_i = c_i_gbox.lower(0);
+   int jmin_c_i = c_i_gbox.lower(1);
+   int jp_c_i = c_i_gbox.numberCells(0);
+   int kmin_c_i = 0;
+   int kp_c_i = 0;
+#if (NDIM == 3)
+   kmin_c_i = c_i_gbox.lower(2);
+   kp_c_i = jp_c_i * c_i_gbox.numberCells(1);
+#endif
+
+   int imin = pbox.lower(0);
+   int imax = pbox.upper(0);
+   int jmin = pbox.lower(1);
+   int jmax = pbox.upper(1);
+   int kmin = 0;
+   int kmax = 0;
+#if (NDIM == 3)
+   kmin = pbox.lower(2);
+   kmax = pbox.upper(2);
+#endif
+
+   std::vector<double> rhs_local(norderp);
+
+   std::vector<double*> ptr_rhs(norderp);
+   for (short i = 0; i < norderp; i++)
+      ptr_rhs[i] = rhs->getPointer(i);
+
+   std::vector<double*> ptr_phi(norderp);
+   for (short i = 0; i < norderp; i++)
+      ptr_phi[i] = phase->getPointer(i);
+
+   // tbox::plog<<"d_norderp_A = "<<d_norderp_A<<std::endl;
+
+   for (int kk = kmin; kk <= kmax; kk++) {
+      for (int jj = jmin; jj <= jmax; jj++) {
+         for (int ii = imin; ii <= imax; ii++) {
+
+            const int idx_rhs = (ii - imin_rhs) + (jj - jmin_rhs) * jp_rhs +
+                                (kk - kmin_rhs) * kp_rhs;
+
+            const int idx_temp = (ii - imin_temp) + (jj - jmin_temp) * jp_temp +
+                                 (kk - kmin_temp) * kp_temp;
+
+            const int idx_pf = (ii - imin_pf) + (jj - jmin_pf) * jp_pf +
+                               (kk - kmin_pf) * kp_pf;
+
+            const int idx_f_i = (ii - imin_f_i) + (jj - jmin_f_i) * jp_f_i +
+                                (kk - kmin_f_i) * kp_f_i;
+
+            const int idx_c_i = (ii - imin_c_i) + (jj - jmin_c_i) * jp_c_i +
+                                (kk - kmin_c_i) * kp_c_i;
+
+            double t = ptr_temp[idx_temp];
+            double fl = ptr_f_l[idx_f_i];
+            double fa = ptr_f_a[idx_f_i];
+            double fb = ptr_f_b[idx_f_i];
+            double cl = ptr_c_l[idx_c_i];
+            double ca = ptr_c_a[idx_c_i];
+            double cb = ptr_c_b[idx_c_i];
+
+            double muL = d_fenergy_strategy->computeMuL(t, cl);
+
+            //
+            // see Moelans, Acta Mat 59 (2011)
+            //
+
+            // driving forces, assuming muL = muA = muB
+            double dfa = (fa - muL * ca);
+            double dfb = (fb - muL * cb);
+            double dfl = (fl - muL * cl);
+            assert(!std::isnan(dfa));
+
+            // interpolation polynomials
+            double hphiA = 0.;
+            for (short i = 0; i < d_norderp_A; i++) {
+               const double phi = std::max(0., ptr_phi[i][idx_pf]);
+               hphiA += phi * phi;
+            }
+            assert(!std::isnan(hphiA));
+
+            double hphiB = 0.;
+            for (short i = d_norderp_A; i < norderp - 1; i++) {
+               const double phi = std::max(0., ptr_phi[i][idx_pf]);
+               hphiB += phi * phi;
+            }
+            assert(!std::isnan(hphiB));
+
+            const double phiL = std::max(0., ptr_phi[norderp - 1][idx_pf]);
+            double hphil = phiL * phiL;
+
+            // std::cout<<"h = "<<hphil<<", "<<hphiA<<", "<<hphiB<<std::endl;
+            const double sum2 = hphil + hphiA + hphiB;
+            assert(sum2 > 0.);
+            const double sum2inv = 1. / sum2;
+
+            hphil *= sum2inv;
+            hphiA *= sum2inv;
+            hphiB *= sum2inv;
+
+            assert(!std::isnan(hphiA));
+            assert(!std::isnan(hphiB));
+
+            // muL ill-defined for hphil very small leads to instabilities
+            // replace driving force by stabilization term in this case
+            // that minimizes ||hphil*cl-hphiA*cA-hphiB*cB-c||
+            // const double epsilon = 1.e-4;
+            //            if (hphil < epsilon) {
+            //const double factor = 1.e4;
+            //double c = ptr_c[idx_pf];
+            // dfl = 0.;
+            //double r = hphil * cl + hphiA * ca + hphiB * cb - c;
+            //          dfl += factor * r * cl;
+            //          dfa += factor * r * ca;
+            //          dfb += factor * r * cb;
+            // std::cout<<"dfb = "<<dfb<<std::endl;
+            //            }
+
+            // solid phase A order parameters
+            for (short i = 0; i < d_norderp_A; i++)
+               rhs_local[i] = 2. * std::max(0., ptr_phi[i][idx_pf]) *
+                              ((1. - hphiA) * dfa - hphil * dfl - hphiB * dfb) *
+                              sum2inv;
+            // solid phase B order parameters
+            for (short i = d_norderp_A; i < norderp - 1; i++)
+               rhs_local[i] = 2. * std::max(0., ptr_phi[i][idx_pf]) *
+                              ((1. - hphiB) * dfb - hphil * dfl - hphiA * dfa) *
+                              sum2inv;
+            // liquid phase order parameter
+            rhs_local[norderp - 1] =
+                2. * std::max(0., ptr_phi[norderp - 1][idx_pf]) *
+                ((1. - hphil) * dfl - hphiA * dfa - hphiB * dfb) * sum2inv;
+            for (short i = 0; i < norderp; i++)
+               assert(!std::isnan(rhs_local[i]));
+
+            for (short i = 0; i < norderp; i++)
+               ptr_rhs[i][idx_rhs] -= rhs_local[i];
+         }
+      }
+   }
+}
diff --git a/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.h b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.h
new file mode 100644
index 00000000..8ea8ea34
--- /dev/null
+++ b/source/MultiOrderBinaryThreePhasesDrivingForceStochioAB.h
@@ -0,0 +1,42 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#ifndef MultiOrderBinaryThreePhasesDrivingForceStochioAB_H
+#define MultiOrderBinaryThreePhasesDrivingForceStochioAB_H
+
+#include "FreeEnergyStrategyBinary.h"
+
+#include "SAMRAI/hier/Patch.h"
+
+class MultiOrderBinaryThreePhasesDrivingForceStochioAB
+{
+ public:
+   MultiOrderBinaryThreePhasesDrivingForceStochioAB(
+       FreeEnergyStrategyBinary* fenergy_strategy, const int conc_id,
+       const int norderp_A);
+
+   void addDrivingForce(SAMRAI::hier::Patch& patch, const int temperature_id,
+                        const int phase_id, const int conc_l_id,
+                        const int conc_a_id, const int conc_b_id,
+                        const int f_l_id, const int f_a_id, const int f_b_id,
+                        const int rhs_id);
+
+ private:
+   /*!
+    * Class to evaluate chemical potential in liquid phase
+    */
+   FreeEnergyStrategyBinary* d_fenergy_strategy;
+
+   const int d_conc_id;
+
+   const int d_norderp_A;
+};
+
+#endif
diff --git a/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.cc b/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.cc
index 85c95811..f85dcb29 100644
--- a/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.cc
+++ b/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.cc
@@ -16,8 +16,6 @@
 ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder::
     ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder(
         const int conc_l_id, const int conc_a_id,
-        const Thermo4PFM::EnergyInterpolationType energy_interp_func_type,
-        const Thermo4PFM::ConcInterpolationType conc_interp_func_type,
         std::shared_ptr<tbox::Database> conc_db)
     : EquilibriumPhaseConcentrationsBinaryMultiOrder(conc_l_id, conc_a_id,
                                                      conc_db)
@@ -33,5 +31,6 @@ ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder::
    double Tref = input_db->getDouble("Tref");
 
    d_fenergy.reset(new Thermo4PFM::ParabolicFreeEnergyFunctionsBinary(
-       Tref, coeffL, coeffA, energy_interp_func_type, conc_interp_func_type));
+       Tref, coeffL, coeffA, Thermo4PFM::EnergyInterpolationType::LINEAR,
+       Thermo4PFM::ConcInterpolationType::LINEAR));
 }
diff --git a/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.h b/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.h
index a429a5e3..edfd4e30 100644
--- a/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.h
+++ b/source/ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder.h
@@ -24,8 +24,6 @@ class ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder
  public:
    ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder(
        const int conc_l_id, const int conc_a_id,
-       const Thermo4PFM::EnergyInterpolationType energy_interp_func_type,
-       const Thermo4PFM::ConcInterpolationType conc_interp_func_type,
        std::shared_ptr<tbox::Database> conc_db);
 
    ~ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder() {}
diff --git a/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.cc b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.cc
new file mode 100644
index 00000000..84af0581
--- /dev/null
+++ b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.cc
@@ -0,0 +1,83 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#include "ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h"
+#include "FuncFort.h"
+#include "ParabolicTools.h"
+
+ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB::
+    ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB(
+        const short norderp_A, const int conc_l_id, const int conc_a_id,
+        const int conc_b_id, const QuatModelParameters& model_parameters,
+        std::shared_ptr<tbox::Database> conc_db)
+    : EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases(
+          norderp_A, conc_l_id, conc_a_id, conc_b_id, model_parameters,
+          conc_db),
+      d_model_parameters(model_parameters)
+{
+}
+
+ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB::
+    ~ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB(){};
+
+int ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB::
+    computePhaseConcentrations(const double temp, double* c, double* hphi,
+                               double* x)
+{
+   assert(!std::isnan(x[0]));
+
+   const double epsilon1 = 1.e-4;
+   const double epsilon2 = 2.e-4;
+
+   // initialize to NaN to trigger error if used when not set
+   double xkks = tbox::IEEE::getSignalingNaN();
+   double xeq = tbox::IEEE::getSignalingNaN();
+
+   const double cA = d_model_parameters.getStochioA();
+   const double cB = d_model_parameters.getStochioB();
+
+   if (hphi[0] >= epsilon1) {
+      // solve explicit KKS problem
+      xkks = (c[0] - hphi[1] * cA - hphi[2] * cB) / hphi[0];
+
+#if 0
+      for (short i = 0; i < 3; i++)
+         std::cerr << hphi[i] << ", ";
+      std::cerr << ", c=" << c[0] << std::endl;
+      std::cerr << "x = " << x[0] << std::endl;
+      std::cerr << "xkks = " << xkks << std::endl;
+      std::cerr << "conc[0] - hphi1 * cA - hphi2 * cB = "
+                << c[0] - hphi[1] * cA - hphi[2] * cB << std::endl;
+#endif
+   }
+
+   if (hphi[0] < epsilon2) {
+      xeq = d_model_parameters.ceq_liquid(temp);
+      // std::cout << "xeq = " << xeq << std::endl;
+      assert(!std::isnan(xeq));
+   }
+
+   if (hphi[0] >= epsilon2) {
+      x[0] = xkks;
+   } else if (hphi[0] <= epsilon1) {
+      x[0] = xeq;
+   } else {
+      // map epsilon1 < hphi < epsilon2 to (0,1)
+      double h = (hphi[0] - epsilon1) / (epsilon2 - epsilon1);
+      // mix equilibrium compositions and KKS solution
+      double f =
+          Thermo4PFM::interp_func(Thermo4PFM::EnergyInterpolationType::PBG, h);
+      x[0] = (1. - f) * xeq + f * xkks;
+   }
+   x[1] = cA;
+   x[2] = cB;
+
+   return 0;
+}
diff --git a/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h
new file mode 100644
index 00000000..bc176310
--- /dev/null
+++ b/source/ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h
@@ -0,0 +1,36 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#ifndef included_ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB
+#define included_ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB
+
+#include "EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases.h"
+#include "QuatModelParameters.h"
+
+class ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB
+    : public EquilibriumPhaseConcentrationsBinaryMultiOrderThreePhases
+{
+ public:
+   ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB(
+       const short norderp_A, const int conc_l_id, const int conc_a_id,
+       const int conc_b_id, const QuatModelParameters& model_parameters,
+       std::shared_ptr<tbox::Database> conc_db);
+
+   ~ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB();
+
+ protected:
+   virtual int computePhaseConcentrations(const double t, double* c,
+                                          double* hphi, double* x);
+
+ private:
+   const QuatModelParameters d_model_parameters;
+};
+
+#endif
diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc
index 18b08472..b6467ce1 100644
--- a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc
+++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.cc
@@ -24,13 +24,12 @@ using namespace SAMRAI;
 
 ParabolicFreeEnergyMultiOrderBinaryThreePhase::
     ParabolicFreeEnergyMultiOrderBinaryThreePhase(
-        std::shared_ptr<tbox::Database> input_db,
-        const Thermo4PFM::ConcInterpolationType conc_interp_func_type,
-        const short norderp_A, MolarVolumeStrategy* mvstrategy,
-        const int conc_l_id, const int conc_a_id, const int conc_b_id)
+        std::shared_ptr<tbox::Database> input_db, const short norderp_A,
+        MolarVolumeStrategy* mvstrategy, const int conc_l_id,
+        const int conc_a_id, const int conc_b_id)
     : FreeEnergyStrategyBinary(Thermo4PFM::EnergyInterpolationType::LINEAR,
-                               conc_interp_func_type, conc_l_id, conc_a_id,
-                               conc_b_id, false),
+                               Thermo4PFM::ConcInterpolationType::LINEAR,
+                               conc_l_id, conc_a_id, conc_b_id, false),
       d_mv_strategy(mvstrategy)
 {
    tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhase..."
diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h
index 79adf021..5f3d4ab6 100644
--- a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h
+++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhase.h
@@ -21,10 +21,9 @@ class ParabolicFreeEnergyMultiOrderBinaryThreePhase
 {
  public:
    ParabolicFreeEnergyMultiOrderBinaryThreePhase(
-       std::shared_ptr<tbox::Database> input_db,
-       const Thermo4PFM::ConcInterpolationType conc_interp_func_type,
-       const short norderp_A, MolarVolumeStrategy* mvstrategy,
-       const int conc_l_id, const int conc_a_id, const int conc_b_id);
+       std::shared_ptr<tbox::Database> input_db, const short norderp_A,
+       MolarVolumeStrategy* mvstrategy, const int conc_l_id,
+       const int conc_a_id, const int conc_b_id);
 
    ~ParabolicFreeEnergyMultiOrderBinaryThreePhase();
 
diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.cc b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.cc
new file mode 100644
index 00000000..78c6c10c
--- /dev/null
+++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.cc
@@ -0,0 +1,219 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#include "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h"
+#include "ParabolicTools.h"
+
+#include "SAMRAI/tbox/InputManager.h"
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/hier/Index.h"
+#include "SAMRAI/math/HierarchyCellDataOpsReal.h"
+
+using namespace SAMRAI;
+
+#include <cassert>
+
+//=======================================================================
+
+ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::
+    ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB(
+        std::shared_ptr<tbox::Database> input_db, const short norderp_A,
+        MolarVolumeStrategy* mvstrategy, const int conc_l_id,
+        const int conc_a_id, const int conc_b_id)
+    : FreeEnergyStrategyBinary(Thermo4PFM::EnergyInterpolationType::LINEAR,
+                               Thermo4PFM::ConcInterpolationType::LINEAR,
+                               conc_l_id, conc_a_id, conc_b_id, false),
+      d_mv_strategy(mvstrategy)
+{
+   tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB..."
+              << std::endl;
+
+   assert(norderp_A > 0);
+   assert(conc_b_id >= 0);
+
+   double coeffL[3][2];
+   readParabolicData(input_db, "Liquid", coeffL);
+
+   double coeffA[3][2];
+   readParabolicData(input_db, "PhaseA", coeffA);
+
+   double coeffB[3][2];
+   readParabolicData(input_db, "PhaseB", coeffB);
+
+   double Tref = input_db->getDouble("Tref");
+
+   d_parabolic_fenergy.reset(
+       new Thermo4PFM::ParabolicFreeEnergyFunctionsBinaryThreePhase(
+           Tref, coeffL, coeffA, coeffB,
+           Thermo4PFM::EnergyInterpolationType::LINEAR,
+           Thermo4PFM::ConcInterpolationType::LINEAR));
+
+   d_multiorder_driving_force.reset(
+       new MultiOrderBinaryThreePhasesDrivingForceStochioAB(this, -1,
+                                                            norderp_A));
+
+   // conversion factor from [J/mol] to [pJ/(mu m)^3]
+   // vm^-1 [mol/m^3] * 10e-18 [m^3/(mu m^3)] * 10e12 [pJ/J]
+   // d_jpmol2pjpmumcube = 1.e-6 / d_vm;
+
+   // R = 8.314472 J · K-1 · mol-1
+   // tbox::plog << "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB:" <<
+   // std::endl; tbox::plog << "Molar volume L =" << vml << std::endl;
+   // tbox::plog << "Molar volume A =" << vma << std::endl;
+   // tbox::plog << "jpmol2pjpmumcube=" << d_jpmol2pjpmumcube << std::endl;
+}
+
+ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::
+    ~ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB(){};
+
+bool ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeCeqT(
+    const double temperature, const Thermo4PFM::PhaseIndex pi0,
+    const Thermo4PFM::PhaseIndex pi1, double* ceq)
+{
+   TBOX_ERROR(
+       "ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeCeqT() "
+       "not "
+       "implemented");
+   return false;
+}
+
+
+//=======================================================================
+
+double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::
+    computeFreeEnergy(const double temperature, double* c_i,
+                      const Thermo4PFM::PhaseIndex pi, const bool gp)
+{
+   assert(d_mv_strategy != nullptr);
+
+   double f = d_parabolic_fenergy->computeFreeEnergy(temperature, c_i, pi, gp);
+   f *= d_mv_strategy->computeInvMolarVolume(temperature, c_i, pi);
+   return f;
+}
+
+//=======================================================================
+
+double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::
+    computeDerivFreeEnergy(const double temperature, double* c_i,
+                           const Thermo4PFM::PhaseIndex pi)
+{
+   double deriv;
+   d_parabolic_fenergy->computeDerivFreeEnergy(temperature, c_i, pi, &deriv);
+   deriv *= d_mv_strategy->computeInvMolarVolume(temperature, c_i, pi);
+   return deriv;
+}
+
+
+void ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::addDrivingForce(
+    const double time, hier::Patch& patch, const int temperature_id,
+    const int phase_id, const int conc_id, const int f_l_id, const int f_a_id,
+    const int f_b_id, const int rhs_id)
+{
+   (void)time;
+   (void)f_b_id;
+
+   d_multiorder_driving_force->addDrivingForce(patch, temperature_id, phase_id,
+                                               d_conc_l_id, d_conc_a_id,
+                                               d_conc_b_id, f_l_id, f_a_id,
+                                               f_b_id, rhs_id);
+}
+
+//=======================================================================
+
+double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeMuL(
+    const double t, const double c0)
+{
+   double c = c0;
+   double mu;
+   d_parabolic_fenergy->computeDerivFreeEnergy(t, &c,
+                                               Thermo4PFM::PhaseIndex::phaseL,
+                                               &mu);
+   return mu *
+          d_mv_strategy->computeInvMolarVolume(t, &c,
+                                               Thermo4PFM::PhaseIndex::phaseL);
+}
+
+//=======================================================================
+
+double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeMuA(
+    const double t, const double c0)
+{
+   double c = c0;
+   double mu;
+   d_parabolic_fenergy->computeDerivFreeEnergy(t, &c,
+                                               Thermo4PFM::PhaseIndex::phaseA,
+                                               &mu);
+   return mu *
+          d_mv_strategy->computeInvMolarVolume(t, &c,
+                                               Thermo4PFM::PhaseIndex::phaseA);
+}
+
+//=======================================================================
+
+double ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::computeMuB(
+    const double t, const double c0)
+{
+   double c = c0;
+   double mu;
+   d_parabolic_fenergy->computeDerivFreeEnergy(t, &c,
+                                               Thermo4PFM::PhaseIndex::phaseB,
+                                               &mu);
+   return mu *
+          d_mv_strategy->computeInvMolarVolume(t, &c,
+                                               Thermo4PFM::PhaseIndex::phaseB);
+}
+
+//=======================================================================
+
+void ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::
+    computeSecondDerivativeEnergyPhaseL(const double temp,
+                                        const std::vector<double>& c_l,
+                                        std::vector<double>& d2fdc2,
+                                        const bool use_internal_units)
+{
+   d_parabolic_fenergy->computeSecondDerivativeFreeEnergy(
+       0., &c_l[0], Thermo4PFM::PhaseIndex::phaseL, &d2fdc2[0]);
+   if (use_internal_units)
+      for (short i = 0; i < 3; i++)
+         d2fdc2[i] *= d_mv_strategy->computeInvMolarVolume(
+             temp, &c_l[0], Thermo4PFM::PhaseIndex::phaseL);
+}
+
+//=======================================================================
+
+void ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::
+    computeSecondDerivativeEnergyPhaseA(const double temp,
+                                        const std::vector<double>& c_a,
+                                        std::vector<double>& d2fdc2,
+                                        const bool use_internal_units)
+{
+   d_parabolic_fenergy->computeSecondDerivativeFreeEnergy(
+       0., &c_a[0], Thermo4PFM::PhaseIndex::phaseA, &d2fdc2[0]);
+   if (use_internal_units)
+      for (short i = 0; i < 3; i++)
+         d2fdc2[i] *= d_mv_strategy->computeInvMolarVolume(
+             temp, &c_a[0], Thermo4PFM::PhaseIndex::phaseA);
+}
+
+//=======================================================================
+
+void ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB::
+    computeSecondDerivativeEnergyPhaseB(const double temp,
+                                        const std::vector<double>& c_b,
+                                        std::vector<double>& d2fdc2,
+                                        const bool use_internal_units)
+{
+   d_parabolic_fenergy->computeSecondDerivativeFreeEnergy(
+       0., &c_b[0], Thermo4PFM::PhaseIndex::phaseB, &d2fdc2[0]);
+   if (use_internal_units)
+      for (short i = 0; i < 3; i++)
+         d2fdc2[i] *= d_mv_strategy->computeInvMolarVolume(
+             temp, &c_b[0], Thermo4PFM::PhaseIndex::phaseB);
+}
diff --git a/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h
new file mode 100644
index 00000000..d9e6b79d
--- /dev/null
+++ b/source/ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB.h
@@ -0,0 +1,73 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#ifndef included_ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB
+#define included_ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB
+
+#include "ParabolicFreeEnergyFunctionsBinaryThreePhase.h"
+#include "FreeEnergyStrategyBinary.h"
+#include "MolarVolumeStrategy.h"
+#include "MultiOrderBinaryThreePhasesDrivingForceStochioAB.h"
+
+class ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB
+    : public FreeEnergyStrategyBinary
+{
+ public:
+   ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB(
+       std::shared_ptr<SAMRAI::tbox::Database> input_db, const short norderp_A,
+       MolarVolumeStrategy* mvstrategy, const int conc_l_id,
+       const int conc_a_id, const int conc_b_id);
+
+   ~ParabolicFreeEnergyMultiOrderBinaryThreePhaseStochioAB();
+
+   void addDrivingForce(const double time, hier::Patch& patch,
+                        const int temperature_id, const int phase_id,
+                        const int conc_id, const int f_l_id, const int f_a_id,
+                        const int f_b_id, const int rhs_id) override;
+
+   void preRunDiagnostics(const double temperature){};
+
+   bool computeCeqT(const double temperature, const Thermo4PFM::PhaseIndex pi0,
+                    const Thermo4PFM::PhaseIndex pi1, double* ceq);
+
+ protected:
+   double computeFreeEnergy(const double temperature, double* c_i,
+                            const Thermo4PFM::PhaseIndex pi,
+                            const bool gp) override;
+
+   double computeDerivFreeEnergy(const double temperature, double* c_i,
+                                 const Thermo4PFM::PhaseIndex pi) override;
+
+ private:
+   MolarVolumeStrategy* d_mv_strategy;
+
+   std::shared_ptr<Thermo4PFM::ParabolicFreeEnergyFunctionsBinaryThreePhase>
+       d_parabolic_fenergy;
+
+   std::shared_ptr<MultiOrderBinaryThreePhasesDrivingForceStochioAB>
+       d_multiorder_driving_force;
+
+   void computeSecondDerivativeEnergyPhaseL(
+       const double temp, const std::vector<double>& c,
+       std::vector<double>& d2fdc2, const bool use_internal_units) override;
+   void computeSecondDerivativeEnergyPhaseA(
+       const double temp, const std::vector<double>& c,
+       std::vector<double>& d2fdc2, const bool use_internal_units) override;
+   void computeSecondDerivativeEnergyPhaseB(const double temp,
+                                            const std::vector<double>& c,
+                                            std::vector<double>& d2fdc2,
+                                            const bool use_internal_units);
+
+   double computeMuL(const double t, const double c);
+   double computeMuA(const double t, const double c);
+   double computeMuB(const double t, const double c);
+};
+
+#endif
diff --git a/source/PhaseConcentrationsStrategyFactory.h b/source/PhaseConcentrationsStrategyFactory.h
index c8bbdbe4..ff527063 100644
--- a/source/PhaseConcentrationsStrategyFactory.h
+++ b/source/PhaseConcentrationsStrategyFactory.h
@@ -20,6 +20,7 @@
 #include "QuadraticEquilibriumPhaseConcentrationsBinaryMultiOrder.h"
 #include "QuadraticEquilibriumThreePhasesTernaryMultiOrder.h"
 #include "ParabolicEquilibriumThreePhasesBinaryMultiOrder.h"
+#include "ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB.h"
 #include "PartitionPhaseConcentrationsStrategy.h"
 #include "PhaseIndependentConcentrationsStrategy.h"
 #include "Database2JSON.h"
@@ -29,14 +30,12 @@
 #include "CALPHADequilibriumPhaseConcentrationsStrategyMultiOrder.h"
 #include "CALPHADequilibriumPhaseConcentrationsStrategyMultiOrderThreePhases.h"
 #include "CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioB.h"
+#include "CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB.h"
 #include "CALPHADFunctions.h"
 
 #include <boost/property_tree/ptree.hpp>
 #include <boost/property_tree/json_parser.hpp>
 
-#include <boost/property_tree/json_parser.hpp>
-#include <boost/property_tree/ptree.hpp>
-
 class PhaseConcentrationsStrategyFactory
 {
  public:
@@ -68,6 +67,11 @@ class PhaseConcentrationsStrategyFactory
 
       std::shared_ptr<PhaseConcentrationsStrategy> phase_conc_strategy;
 
+      const double cA = model_parameters.getStochioA();
+      const double cB = model_parameters.getStochioB();
+      tbox::plog << "cA = " << cA << std::endl;
+      tbox::plog << "cB = " << cB << std::endl;
+
       if (model_parameters.kks_phase_concentration()) {
          tbox::plog << "Phase concentration determined by KKS" << std::endl;
          if (model_parameters.isConcentrationModelCALPHAD()) {
@@ -75,16 +79,21 @@ class PhaseConcentrationsStrategyFactory
             if (ncompositions == 1) {
                tbox::plog << "Binary alloy..." << std::endl;
                const bool subl = Thermo4PFM::checkSublattice(calphad_pt);
-               const double cB = model_parameters.getStochioB();
-               tbox::plog << "cB = " << cB << std::endl;
                if (conc_b_scratch_id >= 0) {
                   tbox::plog << "Three phases..." << std::endl;
                   // three phases
                   if (model_parameters.withMultipleOrderP()) {
                      // multi-order parameters model
                      tbox::plog << "Multi-order parameters..." << std::endl;
-                     if (cB >= 0.) {
-                        tbox::plog << "Stochiometric case..." << std::endl;
+                     if (cB >= 0. && cA >= 0.) {
+                        tbox::plog << "AB stochiometric case..." << std::endl;
+                        phase_conc_strategy.reset(
+                            new CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioAB(
+                                model_parameters.norderpA(), conc_l_scratch_id,
+                                conc_a_scratch_id, conc_b_scratch_id,
+                                model_parameters, conc_db));
+                     } else if (cB >= 0.) {
+                        tbox::plog << "B stochiometric case..." << std::endl;
                         phase_conc_strategy.reset(
                             new CALPHADequilibriumPhaseConcentrationsMultiOrderThreePhasesStochioB(
                                 model_parameters.norderpA(), conc_l_scratch_id,
@@ -213,18 +222,23 @@ class PhaseConcentrationsStrategyFactory
                   if (conc_b_scratch_id >= 0) {
                      tbox::plog << "Parabolic MultiOrder, Three phases..."
                                 << std::endl;
-                     phase_conc_strategy.reset(
-                         new ParabolicEquilibriumThreePhasesBinaryMultiOrder(
-                             model_parameters.norderpA(), conc_l_scratch_id,
-                             conc_a_scratch_id, conc_b_scratch_id,
-                             model_parameters, conc_db));
+                     if (cB >= 0. && cA >= 0.) {
+                        phase_conc_strategy.reset(
+                            new ParabolicEquilibriumThreePhasesBinaryMultiOrderStochioAB(
+                                model_parameters.norderpA(), conc_l_scratch_id,
+                                conc_a_scratch_id, conc_b_scratch_id,
+                                model_parameters, conc_db));
+                     } else {
+                        phase_conc_strategy.reset(
+                            new ParabolicEquilibriumThreePhasesBinaryMultiOrder(
+                                model_parameters.norderpA(), conc_l_scratch_id,
+                                conc_a_scratch_id, conc_b_scratch_id,
+                                model_parameters, conc_db));
+                     }
                   } else {
                      phase_conc_strategy.reset(
                          new ParabolicEquilibriumPhaseConcentrationsBinaryMultiOrder(
-                             conc_l_scratch_id, conc_a_scratch_id,
-                             model_parameters.energy_interp_func_type(),
-                             model_parameters.conc_interp_func_type(),
-                             conc_db));
+                             conc_l_scratch_id, conc_a_scratch_id, conc_db));
                   }
                } else {
                   phase_conc_strategy.reset(
diff --git a/source/QuatIntegrator.cc b/source/QuatIntegrator.cc
index 43b8f43e..fb964d4c 100644
--- a/source/QuatIntegrator.cc
+++ b/source/QuatIntegrator.cc
@@ -3077,7 +3077,8 @@ void QuatIntegrator::setCoefficients(
 
    if (d_with_concentration &&
        d_model_parameters.concentrationModelNeedsPhaseConcentrations()) {
-      d_quat_model->computePhaseConcentrations(hierarchy);
+      d_quat_model->computePhaseConcentrations(time - d_current_time,
+                                               hierarchy);
    }
 
    // mobilities may depend on cl and cs, thus they should be computed after
@@ -3206,7 +3207,8 @@ int QuatIntegrator::evaluateRHSFunction(double time, SundialsAbstractVector* y,
 
             // compute phase concentrations again if they depend on velocity
             // tbox::pout<<"Evaluate c_L, c_S..."<<endl;
-            d_quat_model->computePhaseConcentrations(hierarchy);
+            d_quat_model->computePhaseConcentrations(time - d_current_time,
+                                                     hierarchy);
          }
 
       } while (need_iterate);
diff --git a/source/QuatModel.cc b/source/QuatModel.cc
index d4464fde..a7f2114e 100644
--- a/source/QuatModel.cc
+++ b/source/QuatModel.cc
@@ -347,8 +347,8 @@ void QuatModel::initializeRHSandEnergyStrategies(
    d_free_energy_strategy =
        FreeEnergyStrategyFactory::create(d_model_parameters, d_ncompositions,
                                          d_conc_l_id, d_conc_a_id, d_conc_b_id,
-                                         d_mvstrategy, d_meltingT_strategy,
-                                         Tref, d_conc_db);
+                                         d_conc_scratch_id, d_mvstrategy,
+                                         d_meltingT_strategy, Tref, d_conc_db);
 
    if (d_model_parameters.with_concentration()) {
       if (d_model_parameters.concentrationModelNeedsPhaseConcentrations()) {
@@ -734,6 +734,7 @@ void QuatModel::Initialize(std::shared_ptr<tbox::MemoryDatabase>& input_db,
    const double cLref = d_model_parameters.concL_ref();
    if (cLref >= 0.) {
       tbox::plog << "Set reference cL to " << cLref << std::endl;
+      assert(cLref <= 1.);
       if (d_conc_l_ref_id >= 0)
          cellops.setToScalar(d_conc_l_ref_id, cLref, false);
       cellops.setToScalar(d_conc_l_id, cLref, false);
@@ -1213,7 +1214,8 @@ void QuatModel::registerPhaseConcentrationVariables()
       assert(d_conc_b_var);
    }
 
-   if (d_model_parameters.isConcentrationModelCALPHAD()) {
+   if (d_model_parameters.isConcentrationModelCALPHAD() ||
+       d_model_parameters.getStochioA() >= 0.) {
       d_conc_l_ref_var.reset(
           new pdat::CellVariable<double>(tbox::Dimension(NDIM), "conc_l_ref",
                                          d_ncompositions));
@@ -2323,7 +2325,8 @@ double QuatModel::Advance(void)
    bool update_solution = false;
 
    if (d_model_parameters.with_concentration() &&
-       d_model_parameters.isConcentrationModelCALPHAD())
+       (d_model_parameters.isConcentrationModelCALPHAD() ||
+        d_model_parameters.getStochioA() >= 0.))
       resetRefPhaseConcentrations();
 
    if (d_model_parameters.with_orientation()) {
@@ -4937,12 +4940,9 @@ void QuatModel::evaluateEnergy(
          if (d_model_parameters.needGhosts4PartitionCoeff())
             fillPartitionCoeffGhosts();
       }
-      if (d_model_parameters.concentrationModelNeedsPhaseConcentrations()) {
-         assert(d_phase_conc_strategy != nullptr);
-         d_phase_conc_strategy->computePhaseConcentrations(
-             hierarchy, d_temperature_scratch_id, d_phase_scratch_id,
-             d_conc_scratch_id);
-      }
+      // if (d_model_parameters.concentrationModelNeedsPhaseConcentrations()) {
+      //   computePhaseConcentrations(0., hierarchy);
+      //}
    }
 
    if (d_free_energy_strategy) {
@@ -4984,7 +4984,7 @@ void QuatModel::evaluateEnergy(
 //=======================================================================
 
 void QuatModel::computePhaseConcentrations(
-    const std::shared_ptr<hier::PatchHierarchy> hierarchy)
+    const double dt, const std::shared_ptr<hier::PatchHierarchy> hierarchy)
 {
    assert(d_phase_conc_strategy != nullptr);
 
@@ -5336,6 +5336,8 @@ void QuatModel::setRefPhaseConcentrationsToEquilibrium(const double* const ceq)
 {
    assert(d_conc_l_ref_id >= 0);
    assert(d_conc_a_ref_id >= 0);
+   assert(ceq[0] >= 0.);
+   assert(ceq[0] <= 1.);
 
    tbox::pout << "QuatModel::setRefPhaseConcentrationsToEquilibrium(ceq)"
               << std::endl;
diff --git a/source/QuatModel.h b/source/QuatModel.h
index 806619f9..888b51cd 100644
--- a/source/QuatModel.h
+++ b/source/QuatModel.h
@@ -225,7 +225,7 @@ class QuatModel : public PFModel
        std::shared_ptr<hier::PatchHierarchy> hierarchy);
 
    void computePhaseConcentrations(
-       std::shared_ptr<hier::PatchHierarchy> hierarchy);
+       const double time, std::shared_ptr<hier::PatchHierarchy> hierarchy);
 
    void evaluateEnergy(const std::shared_ptr<hier::PatchHierarchy> hierarchy,
                        const double time, double& total_energy,
diff --git a/source/QuatModelParameters.cc b/source/QuatModelParameters.cc
index 5006ab49..d93fc97b 100644
--- a/source/QuatModelParameters.cc
+++ b/source/QuatModelParameters.cc
@@ -90,6 +90,7 @@ QuatModelParameters::QuatModelParameters() : d_moving_frame_velocity(def_val)
    d_liquidus_slope = def_val;
    d_average_concentration = def_val;
    d_stochio_cB = -1.;
+   d_stochio_cA = -1.;
 
    for (short i = 0; i < 3; i++) {
       d_ceq0_solidA[i] = def_val;
@@ -433,6 +434,7 @@ void QuatModelParameters::readConcDB(std::shared_ptr<tbox::Database> conc_db)
        conc_db->getBoolWithDefault("init_phase_conc_eq", true);
 
    d_stochio_cB = conc_db->getDoubleWithDefault("stochio_cB", -1.);
+   d_stochio_cA = conc_db->getDoubleWithDefault("stochio_cA", -1.);
 
    readEquilibriumCompositions(conc_db);
 }
diff --git a/source/QuatModelParameters.h b/source/QuatModelParameters.h
index 31bbc23c..b7798414 100644
--- a/source/QuatModelParameters.h
+++ b/source/QuatModelParameters.h
@@ -141,7 +141,10 @@ class QuatModelParameters
    {
       return d_initc_in_phase[d_ncompositions + index];
    }
+
    double getStochioB() const { return d_stochio_cB; }
+   double getStochioA() const { return d_stochio_cA; }
+
    double meltingT() const { return d_meltingT; }
    double interfaceMobility() const { return d_interface_mobility; }
    double rescale_factorT() const { return d_rescale_factorT; }
@@ -626,6 +629,7 @@ class QuatModelParameters
 
    // stochio composition of phase B (if set to value >= 0)
    double d_stochio_cB;
+   double d_stochio_cA;
 
    // free energy parameters:
    // f(phi) = d_phase_well_scale * g(phi)
diff --git a/source/fortran/ConcFort.h b/source/fortran/ConcFort.h
index 427022e9..51d23e8e 100644
--- a/source/fortran/ConcFort.h
+++ b/source/fortran/ConcFort.h
@@ -183,6 +183,17 @@ void ADD_FLUX_4TH(const int& ifirst0, const int& ilast0, const int& ifirst1,
 #endif
                   const int& ngflux, const int* physb);
 
+void ZERO_FLUX_PHYSB(const int& ifirst0, const int& ilast0, const int& ifirst1,
+                     const int& ilast1,
+#if (NDIM == 3)
+                     const int& ifirst2, const int& ilast2,
+#endif
+                     const double* flux0, const double* flux1,
+#if (NDIM == 3)
+                     const double* flux2,
+#endif
+                     const int& ngflux, const int* physb);
+
 void CONCENTRATION_FLUX_SPINODAL(
     const int& ifirst0, const int& ilast0, const int& ifirst1,
     const int& ilast1,