diff --git a/samples/matlab/lithium_ion_battery.m b/samples/matlab/lithium_ion_battery.m
index 200c67105..0113cf485 100644
--- a/samples/matlab/lithium_ion_battery.m
+++ b/samples/matlab/lithium_ion_battery.m
@@ -48,12 +48,12 @@ S_ca = 1.1167; % [m^2] Cathode total active material surface area
 S_an = 0.7824; % [m^2] Anode total active material surface area
 
 % Import all Cantera phases
-anode = importThermoPhase(inputCTI, 'anode');
-cathode = importThermoPhase(inputCTI,'cathode');
-elde = importThermoPhase(inputCTI,'electron');
-elyt = importThermoPhase(inputCTI,'electrolyte');
-anode_interface = importEdge(inputCTI, 'edge_anode_electrolyte', anode, elde, elyt);
-cathode_interface = importEdge(inputCTI, 'edge_cathode_electrolyte', cathode, elde, elyt);
+anode = Solution(inputCTI, 'anode');
+cathode = Solution(inputCTI, 'cathode');
+elde = Solution(inputCTI, 'electron');
+elyt = Solution(inputCTI, 'electrolyte');
+anode_interface = Interface(inputCTI, 'edge_anode_electrolyte', anode, elde, elyt);
+cathode_interface = Interface(inputCTI, 'edge_cathode_electrolyte', cathode, elde, elyt);
 
 % Set the temperatures and pressures of all phases
 phases = [anode elde elyt cathode];
@@ -91,13 +91,6 @@ set(gca,'fontsize',14)
 %--------------------------------------------------------------------------
 % Helper functions
 
-% This function returns the ThermoPhase class instance from CTI file
-function phase = importThermoPhase(inputCTI, name)
-    doc = XML_Node('doc', inputCTI);
-    node = findByID(doc, name);
-    phase = ThermoPhase(node);
-end
-
 % This function returns the Cantera calculated anode current (in A)
 function anCurr = anode_curr(phi_s,phi_l,X_Li_an,anode,elde,elyt,anode_interface,S_an)
 
diff --git a/src/thermo/IdealSolidSolnPhase.cpp b/src/thermo/IdealSolidSolnPhase.cpp
index 7041397d6..f82d4d472 100644
--- a/src/thermo/IdealSolidSolnPhase.cpp
+++ b/src/thermo/IdealSolidSolnPhase.cpp
@@ -98,7 +98,7 @@ void IdealSolidSolnPhase::setDensity(const doublereal rho)
     // Unless the input density is exactly equal to the density calculated and
     // stored in the State object, we throw an exception. This is because the
     // density is NOT an independent variable.
-    if (rho != density()) {
+    if (std::abs(rho/density() - 1.0) > 1e-15) {
         throw CanteraError("IdealSolidSolnPhase::setDensity",
                            "Density is not an independent variable");
     }
@@ -365,6 +365,7 @@ bool IdealSolidSolnPhase::addSpecies(shared_ptr<Species> spec)
             throw CanteraError("IdealSolidSolnPhase::addSpecies",
                 "Molar volume not specified for species '{}'", spec->name);
         }
+        calcDensity();
     }
     return added;
 }