diff --git a/src/specmicp/adimensional/adimensional_system_solution_extractor.cpp b/src/specmicp/adimensional/adimensional_system_solution_extractor.cpp
index 40ccb12..be6cfc2 100644
--- a/src/specmicp/adimensional/adimensional_system_solution_extractor.cpp
+++ b/src/specmicp/adimensional/adimensional_system_solution_extractor.cpp
@@ -1,75 +1,80 @@
 #include "adimensional_system_solution_extractor.hpp"
 #include "physics/laws.hpp"
 
 namespace specmicp {
 
+scalar_t AdimensionalSystemSolutionExtractor::density_water()
+{
+    return laws::density_water(units::celsius(25.0), length_unit(), mass_unit());
+}
+
 scalar_t AdimensionalSystemSolutionExtractor::mass_concentration_water()
 {
-    return laws::density_water(units::celsius(25.0), length_unit(), mass_unit())*total_saturation_water();
+    return density_water()*total_saturation_water();
 }
 
 scalar_t AdimensionalSystemSolutionExtractor::pH()
 {
     // find species responsible for pH
     specmicp::database::Database db_handler(m_data);
     index_t id = db_handler.component_label_to_id("HO[-]");
     if (id != no_species)
     {
         return 14+log_activity_component(id);
     }
     else
     {
         id = db_handler.component_label_to_id("H[+]");
         if (id != no_species)
             return -log_activity_component(id);
         throw std::runtime_error("No component corresponding to the dissociation of water !");
     }
 }
 
 scalar_t AdimensionalSystemSolutionExtractor::total_saturation_minerals()
 {
     return m_solution.main_variables.segment(m_data->nb_component, m_data->nb_mineral).sum();
 }
 
 scalar_t AdimensionalSystemSolutionExtractor::mole_concentration_mineral(index_t mineral)
 {
     return total_saturation_mineral(mineral)/m_data->molar_volume_mineral(mineral, length_unit());
 }
 
 scalar_t AdimensionalSystemSolutionExtractor::mass_concentration_mineral(index_t mineral)
 {
     return mole_concentration_mineral(mineral)*m_data->molar_mass_mineral(mineral, mass_unit());
 }
 
 //! \brief Return the total aqueous concentration
 scalar_t AdimensionalSystemSolutionExtractor::total_aqueous_concentration(index_t component)
 {
     scalar_t conc = molality_component(component);
     for (index_t aqueous: m_data->range_aqueous())
     {
         if (m_data->nu_aqueous(aqueous, component) != 0.0)
         {
             conc += m_data->nu_aqueous(aqueous, component)*molality_aqueous(aqueous);
         }
     }
     return conc;
 }
 
 //! \brief Return the total solid concentration
 scalar_t AdimensionalSystemSolutionExtractor::total_solid_concentration(index_t component)
 {
     scalar_t conc= 0;
     for (index_t mineral: m_data->range_mineral())
     {
         if (m_data->nu_mineral(mineral, component) != 0.0)
         {
             conc +=   m_data->nu_mineral(mineral, component)
                     * total_saturation_mineral(mineral)
                     / m_data->molar_volume_mineral(mineral, length_unit());
         }
     }
     return conc;
 }
 
 
 } // end namespace specmicp
diff --git a/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp b/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp
index a9c72f5..7d5182e 100644
--- a/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp
+++ b/src/specmicp/adimensional/adimensional_system_solution_extractor.hpp
@@ -1,194 +1,198 @@
 #ifndef SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONEXTRACTOR_HPP
 #define SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONEXTRACTOR_HPP
 
 #include "common.hpp"
 #include "database.hpp"
 #include "adimensional_system_solution.hpp"
 #include "physics/units.hpp"
 
 //! \file adimensional_system_solution_extractor.hpp Obtain data from AdimensionalSystemSolution
 
 namespace specmicp {
 
 //! \brief Obtain physical variables from AdimensionalSystemSolution
 class AdimensionalSystemSolutionExtractor: public units::UnitBaseClass
 {
 public:
     // public members
     // ###############
 
     AdimensionalSystemSolutionExtractor(
             const AdimensionalSystemSolution& solution,
             RawDatabasePtr thedatabase,
             units::UnitsSet units_set):
         units::UnitBaseClass(units_set),
         m_solution(solution),
         m_data(thedatabase)
     {}
 
     // Primary variables
     // =================
 
     // water
     // ------
 
     //! \brief Return the total saturation of water
     scalar_t total_saturation_water() {
         return m_solution.main_variables(id_water());}
 
+    scalar_t density_water();
+
     // component
     // ---------
 
     //! \brief Return the log_10 of the molality of 'component'
     scalar_t log_molality_component(index_t component) {
         return m_solution.main_variables(id_component(component));}
 
     //! \brief Return the molality of 'component'
     scalar_t molality_component(index_t component) {
+        if (component == 0)
+            return 1.0/m_data->molar_mass_basis(component, mass_unit());
         return pow10(log_molality_component(component));
     }
 
     // solid phases
     // ------------
 
     //! \brief Return the total saturation of 'mineral'
     scalar_t total_saturation_mineral(index_t mineral) {
         return m_solution.main_variables(id_mineral(mineral));
     }
 
     // Secondary variables
     // ===================
 
     //! \brief Return the ionic strength of the solution
     scalar_t ionic_strength() {return m_solution.ionic_strength;}
 
     // component
     // ---------
 
     //! \brief Return the log_10 of the activity coefficient of 'component'
     scalar_t log_activity_coefficient_component(index_t component) {
         return m_solution.log_gamma(id_component_gamma(component));
     }
     //! \brief Return the activity coefficient of 'component'
     scalar_t activity_coefficient_component(index_t component) {
         return pow10(log_activity_coefficient_component(component));
     }
     //! \brief Return the log_10 of the activity of 'component
     scalar_t log_activity_component(index_t component) {
         return log_molality_component(component) + log_activity_coefficient_component(component);
     }
     //! \brief Return the activity of 'component'
     scalar_t activity_component(index_t component) {
         return pow10(log_activity_component(component));
     }
 
     // aqueous species
     // ---------------
     //! \brief Return the molality of secondary specis 'aqueous'
     scalar_t molality_aqueous(index_t aqueous) {
         return m_solution.secondary_molalities(id_aqueous(aqueous));}
     //! \brief Return the log10 of the activity ocefficient of secondary species 'aqueous'
     scalar_t log_activity_coefficient_aqueous(index_t aqueous) {
         return m_solution.log_gamma(id_aqueous_gamma(aqueous));
     }
     //! \brief Return the activity coefficient of secondary species 'aqueous'
     scalar_t activity_coefficient_aqueous(index_t aqueous) {
         return pow10(log_activity_coefficient_aqueous(aqueous));
     }
     //! \brief Return the activity of secondary species 'aqueous'
     scalar_t activity_aqueous(index_t aqueous) {
         return activity_coefficient_aqueous(aqueous)*molality_aqueous(aqueous);
     }
 
     // gas fugacity
     // -------------
     //! \brief Return fugacity for 'gas'
     scalar_t gas_fugacity(index_t gas) {
         return m_solution.gas_fugacities(gas);
     }
 
     // Tertiary variables
     // ==================
 
     // Water
     // -----
 
     //! \brief Return the saturation of water
     scalar_t saturation_water() {return total_saturation_water()*porosity();}
 
     //! \brief Return the mass concentration of water
     scalar_t mass_concentration_water();
 
     //! \brief Return the pH of the solution
     scalar_t pH();
 
     // Component
     // ---------
 
     //! \brief Return the total aqueous concentration
     scalar_t total_aqueous_concentration(index_t component);
     //! \brief Return the total solid concentration
     scalar_t total_solid_concentration(index_t component);
 
     // Gas phase
     // ---------
 
     //! \brief Return the saturation of the gas phase
     scalar_t saturation_gas_phase() {return 1-saturation_water();}
 
     //! \brief Return the total saturation of the gas phase
     scalar_t total_saturation_gas_phase() {return porosity() - saturation_water();}
 
     // Component
     // ---------
 
     // Solid phases
     // ------------
 
     //! \brief Return the concentration of 'mineral'
     scalar_t mole_concentration_mineral(index_t mineral);
 
     //! \brief Return the mass concentration of 'mineral'
     scalar_t mass_concentration_mineral(index_t mineral);
 
     //! \brief Return the total saturation of all minerals
     scalar_t total_saturation_minerals();
 
     //! \brief Return the porosity
     scalar_t porosity() {return 1 - total_saturation_minerals();}
 
 private:
     // private members
     // ###############
 
     index_t id_water() {return 0;}
     index_t id_component(index_t component) {
         specmicp_assert(component < m_data->nb_component);
         return component;
     }
     index_t id_mineral(index_t mineral) {
         specmicp_assert(mineral < m_data->nb_mineral);
         return m_data->nb_component+mineral;}
     index_t id_component_gamma(index_t component) {
         specmicp_assert(component < m_data->nb_component);
         return component;
     }
     index_t id_aqueous(index_t aqueous) {
         specmicp_assert(aqueous < m_data->nb_aqueous);
         return aqueous;
     }
     index_t id_aqueous_gamma(index_t aqueous) {
         specmicp_assert(aqueous < m_data->nb_aqueous);
         return m_data->nb_component+aqueous;
     }
 
 private:
     // private attributes
     // ##################
-    AdimensionalSystemSolution m_solution;
+    const AdimensionalSystemSolution& m_solution;
     RawDatabasePtr m_data;
 };
 
 } // end namespace specmicp
 
 #endif // SPECMICP_SPECMICP_ADIMENSIONALSYSTEMSOLUTIONEXTRACTOR_HPP