diff --git a/src/main/java/neqsim/thermo/phase/PhaseEos.java b/src/main/java/neqsim/thermo/phase/PhaseEos.java
index 642cc7ade..3bfac0e23 100644
--- a/src/main/java/neqsim/thermo/phase/PhaseEos.java
+++ b/src/main/java/neqsim/thermo/phase/PhaseEos.java
@@ -193,10 +193,16 @@ public void setMixingRule(MixingRuleTypeInterface mr) {
   /** {@inheritDoc} */
   @Override
   public void setMixingRuleGEModel(String name) {
-    if (mixRule != null) {
+    if (mixRule == null) {
+      // do nothing or initialize?
+      logger.debug("mixRule is null");
+    } else {
       mixRule.setMixingRuleGEModel(name);
     }
-    if (mixSelect != null) {
+    if (mixSelect == null) {
+      // do nothing or initialize?
+      logger.debug("mixSelect is null");
+    } else {
       mixSelect.setMixingRuleGEModel(name);
     }
   }
@@ -212,7 +218,9 @@ public void resetMixingRule(MixingRuleTypeInterface mr) {
     if (mr == null) {
       mixRule = null;
     } else {
-      mixRule = mixSelect.getMixingRule(mr.getValue(), this);
+      mixRule = mixSelect.resetMixingRule(mr.getValue(), this);
+      // TODO: verify if should resetMixingRule or getMixingRule
+      // mixRule = mixSelect.getMixingRule(mr.getValue(), this);
     }
   }
 
diff --git a/src/main/java/neqsim/thermo/system/SystemThermo.java b/src/main/java/neqsim/thermo/system/SystemThermo.java
index 4a7a40f01..08a73df42 100644
--- a/src/main/java/neqsim/thermo/system/SystemThermo.java
+++ b/src/main/java/neqsim/thermo/system/SystemThermo.java
@@ -259,6 +259,7 @@ public void addComponent(int index, double moles) {
       }
     }
     setTotalNumberOfMoles(getTotalNumberOfMoles() + moles);
+    // TODO: isInitialized = false;
   }
 
   /** {@inheritDoc} */
@@ -280,6 +281,7 @@ public void addComponent(int index, double moles, int phaseNumber) {
     }
 
     setTotalNumberOfMoles(getTotalNumberOfMoles() + moles);
+    // TODO: isInitialized = false;
   }
 
   /** {@inheritDoc} */
@@ -334,6 +336,7 @@ public void addComponent(String componentName, double moles) {
       }
     }
     setTotalNumberOfMoles(getTotalNumberOfMoles() + moles);
+    // TODO: isInitialized = false;
   }
 
   /** {@inheritDoc} */
@@ -359,6 +362,7 @@ public void addComponent(String componentName, double moles, double TC, double P
       componentNames.remove("default");
       componentNames.add(componentName);
     }
+    // TODO: isInitialized = false;
   }
 
   /** {@inheritDoc} */
@@ -398,6 +402,7 @@ public void addComponent(String componentName, double moles, int phaseNumber) {
       getPhase(i).setAttractiveTerm(attractiveTermNumber);
     }
     numberOfComponents++;
+    // TODO: isInitialized = false;
   }
 
   /** {@inheritDoc} */
@@ -428,6 +433,7 @@ public void addComponent(String componentName, double value, String unitName) {
     double SIval = unit.getSIvalue();
     // System.out.println("number of moles " + SIval);
     this.addComponent(componentName, SIval);
+    // TODO: isInitialized = false;
   }
 
   /** {@inheritDoc} */
@@ -458,6 +464,7 @@ public void addComponent(String componentName, double value, String name, int ph
     double SIval = unit.getSIvalue();
     // System.out.println("number of moles " + SIval);
     this.addComponent(componentName, SIval, phaseNum);
+    // TODO: isInitialized = false;
   }
 
   /** {@inheritDoc} */
@@ -863,6 +870,7 @@ public void addTBPfraction(String componentName, double numberOfMoles, double mo
       refSystem.setTemperature(273.15 + 15.0);
       refSystem.setPressure(ThermodynamicConstantsInterface.referencePressure);
       refSystem.addComponent("default", 1.0, 273.15, 50.0, 0.1);
+      refSystem.setMixingRule(1);
       refSystem.init(0);
       refSystem.setNumberOfPhases(1);
       refSystem.setPhaseType(0, PhaseType.LIQUID);
@@ -1697,6 +1705,7 @@ public double getAntoineVaporPressure(double temp) {
   /** {@inheritDoc} */
   @Override
   public final double getBeta() {
+    // TODO: verify, actually returning the heaviest?
     return beta[0];
   }
 
@@ -2507,6 +2516,12 @@ public double[] getMolarRate() {
   /** {@inheritDoc} */
   @Override
   public double getMolarVolume() {
+    if (!this.isInitialized) {
+      this.init(0);
+    }
+    if (!isBetaValid()) {
+      logger.warn("getMolarVolume", "Calculation is wrong, as beta is not valid. Perform flash");
+    }
     double volume = 0;
     for (int i = 0; i < numberOfPhases; i++) {
       volume += beta[phaseIndex[i]] * getPhase(i).getMolarVolume();
@@ -2517,6 +2532,12 @@ public double getMolarVolume() {
   /** {@inheritDoc} */
   @Override
   public double getMolarVolume(String unit) {
+    if (!this.isInitialized) {
+      this.init(0);
+    }
+    if (!isBetaValid()) {
+      logger.warn("getMolarVolume", "Calculation is wrong, as beta is not valid. Perform flash");
+    }
     double volume = 0;
     for (int i = 0; i < numberOfPhases; i++) {
       volume += beta[phaseIndex[i]] * getPhase(i).getMolarVolume(unit);
@@ -2908,6 +2929,16 @@ public final double getSumBeta() {
     return sum;
   }
 
+  /**
+   * Verify if sum of beta is 1. Used to check if System needs to be flashed.
+   *
+   * @return True if the sum of beta is close to 1.
+   */
+  public boolean isBetaValid() {
+    return this.getSumBeta() > 1.0 - ThermodynamicModelSettings.phaseFractionMinimumLimit
+        && this.getSumBeta() < 1.0 + ThermodynamicModelSettings.phaseFractionMinimumLimit;
+  }
+
   /** {@inheritDoc} */
   @Override
   public final double getTC() {
@@ -3307,7 +3338,7 @@ public final void initBeta() {
     if (!isInitialized
         && this.getSumBeta() < 1.0 - ThermodynamicModelSettings.phaseFractionMinimumLimit
         || this.getSumBeta() > 1.0 + ThermodynamicModelSettings.phaseFractionMinimumLimit) {
-      logger.warn("SystemThermo:initBeta - Sum of beta does not equal 1.0 ");
+      logger.warn("SystemThermo:initBeta - Sum of beta does not equal 1.0. " + beta);
     }
   }
 
@@ -3672,7 +3703,6 @@ public void normalizeBeta() {
   @Override
   public void orderByDensity() {
     boolean change = false;
-    // int count = 0;
 
     for (int i = 0; i < getNumberOfPhases(); i++) {
       if (getPhase(i).getPhysicalProperties() == null) {
@@ -3683,9 +3713,9 @@ public void orderByDensity() {
 
     do {
       change = false;
-      // count++;
       for (int i = 1; i < getNumberOfPhases(); i++) {
         if (i == 4) {
+          // Do not sort phase 5 and 6
           break;
         }
 
@@ -3694,7 +3724,7 @@ public void orderByDensity() {
             getPhase(i).initPhysicalProperties(PhysicalPropertyType.MASS_DENSITY);
           }
         } catch (Exception ex) {
-          logger.error(ex.getMessage(), ex);
+          logger.error(ex.getMessage());
         }
         if (getPhase(i).getPhysicalProperties().calcDensity() < getPhase(i - 1)
             .getPhysicalProperties().calcDensity()) {
@@ -4002,6 +4032,7 @@ public void reset() {
       addComponent(getPhase(0).getComponent(i).getComponentName(),
           -getPhase(0).getComponent(i).getNumberOfmoles());
     }
+    // TODO: isInitialized = false;
   }
 
   /** {@inheritDoc} */
@@ -4248,10 +4279,13 @@ public void setAttractiveTerm(int i) {
   /** {@inheritDoc} */
   @Override
   public final void setBeta(double b) {
+    // TODO: if number of phases > 2, should fail
     if (b < 0) {
+      logger.warn("setBeta - Tried to set beta < 0: " + beta);
       b = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
     }
     if (b > 1) {
+      logger.warn("setBeta - Tried to set beta > 1: " + beta);
       b = 1.0 - neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
     }
     beta[0] = b;
@@ -4262,9 +4296,11 @@ public final void setBeta(double b) {
   @Override
   public final void setBeta(int phaseNum, double b) {
     if (b < 0) {
+      logger.warn("setBeta - Tried to set beta < 0: " + beta);
       b = neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
     }
     if (b > 1) {
+      logger.warn("setBeta - Tried to set beta > 1: " + beta);
       b = 1.0 - neqsim.thermo.ThermodynamicModelSettings.phaseFractionMinimumLimit;
     }
     beta[phaseIndex[phaseNum]] = b;