diff --git a/.gitmodules b/.gitmodules
index 272df3f35..f61e8946a 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -10,4 +10,4 @@
 [submodule "extern/Microphysics"]
 	path = extern/Microphysics
 	url = https://github.com/psharda/Microphysics
-	branch = del-eos-comp
+	branch = development
diff --git a/extern/Microphysics b/extern/Microphysics
index 5b4b3feee..ea4204146 160000
--- a/extern/Microphysics
+++ b/extern/Microphysics
@@ -1 +1 @@
-Subproject commit 5b4b3feee521f460d46401e80f34200846a1545f
+Subproject commit ea4204146eb2269605764e432a022cb1ae090398
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a6b1db367..18ef20512 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -95,7 +95,9 @@ link_libraries(hdf5::hdf5)
 include(CTest)
 
 #create an object library for files that should be compiled with all test problems
-set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/CloudyCooling.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/GrackleDataReader.cpp")
+set (QuokkaObjSources "${CMAKE_CURRENT_SOURCE_DIR}/main.cpp" "${CMAKE_CURRENT_SOURCE_DIR}/CloudyCooling.cpp" 
+                      "${CMAKE_CURRENT_SOURCE_DIR}/GrackleDataReader.cpp"
+                      "${gamma_law_sources}")
 #we don't use it anymore because it gives issues on Cray systems
 #add_library(QuokkaObj OBJECT ${QuokkaObjSources} ${gamma_law_sources})
 #if(AMReX_GPU_BACKEND MATCHES "CUDA")
diff --git a/src/CloudyCooling.hpp b/src/CloudyCooling.hpp
index 2e5dfe50c..ae5aca447 100644
--- a/src/CloudyCooling.hpp
+++ b/src/CloudyCooling.hpp
@@ -81,7 +81,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto cloudy_cooling_function(Real const
 {
 	// interpolate cooling rates from Cloudy tables
 	const Real rhoH = rho * cloudy_H_mass_fraction; // mass density of H species
-	const Real nH = rhoH / quokka::hydrogen_mass_cgs;
+	const Real nH = rhoH / (C::m_p + C::m_e);
 	const Real log_nH = std::log10(nH);
 	const Real log_T = std::log10(T);
 
@@ -102,7 +102,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto cloudy_cooling_function(Real const
 
 	// compute electron density
 	// N.B. it is absolutely critical to include the metal contribution here!
-	double n_e = (rho / quokka::hydrogen_mass_cgs) * (1.0 - mu * (X + Y / 4. + Z / mean_metals_A)) / (mu - (electron_mass_cgs / quokka::hydrogen_mass_cgs));
+	double n_e = (rho / (C::m_p + C::m_e)) * (1.0 - mu * (X + Y / 4. + Z / mean_metals_A)) / (mu - (electron_mass_cgs / (C::m_p + C::m_e)));
 	// the approximation for the metals contribution to e- fails at high densities (~1e3 or higher)
 	n_e = std::max(n_e, 1.0e-4 * nH);
 
@@ -118,7 +118,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto cloudy_cooling_function(Real const
 	// Compton term (CMB photons)
 	// [e.g., Hirata 2018: doi:10.1093/mnras/stx2854]
 	constexpr double Gamma_C = (8. * sigma_T * E_cmb) / (3. * electron_mass_cgs * c_light_cgs_);
-	constexpr double C_n = Gamma_C * quokka::boltzmann_constant_cgs / (5. / 3. - 1.0);
+	constexpr double C_n = Gamma_C * C::k_B / (5. / 3. - 1.0);
 	const double compton_CMB = -C_n * (T - T_cmb) * n_e;
 	Edot += compton_CMB;
 
@@ -129,13 +129,13 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto ComputeEgasFromTgas(double rho, do
 {
 	// convert Egas (internal gas energy) to temperature
 	const Real rhoH = rho * cloudy_H_mass_fraction;
-	const Real nH = rhoH / quokka::hydrogen_mass_cgs;
+	const Real nH = rhoH / (C::m_p + C::m_e);
 
 	// compute mu from mu(T) table
 	const Real mu = interpolate2d(std::log10(nH), std::log10(Tgas), tables.log_nH, tables.log_Tgas, tables.meanMolWeight);
 
 	// compute thermal gas energy
-	const Real Egas = (rho / (quokka::hydrogen_mass_cgs * mu)) * quokka::boltzmann_constant_cgs * Tgas / (gamma - 1.);
+	const Real Egas = (rho / ((C::m_p + C::m_e) * mu)) * C::k_B * Tgas / (gamma - 1.);
 	return Egas;
 }
 
@@ -157,12 +157,12 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto ComputeTgasFromEgas(double rho, do
 
 	// solve for temperature given Eint (with fixed adiabatic index gamma)
 	const Real rhoH = rho * cloudy_H_mass_fraction;
-	const Real nH = rhoH / quokka::hydrogen_mass_cgs;
+	const Real nH = rhoH / (C::m_p + C::m_e);
 	const Real log_nH = std::log10(nH);
 
 	// mean molecular weight (in Grackle tables) is defined w/r/t
 	// hydrogen_mass_cgs_
-	const Real C = (gamma - 1.) * Egas / (quokka::boltzmann_constant_cgs * (rho / quokka::hydrogen_mass_cgs));
+	const Real C = (gamma - 1.) * Egas / (C::k_B * (rho / (C::m_p + C::m_e)));
 
 	// solve for mu(T)*C == T.
 	// (Grackle does this with a fixed-point iteration. We use a more robust
diff --git a/src/Cooling/test_cooling.cpp b/src/Cooling/test_cooling.cpp
index e602c79cf..e88164691 100644
--- a/src/Cooling/test_cooling.cpp
+++ b/src/Cooling/test_cooling.cpp
@@ -30,13 +30,13 @@ using amrex::Real;
 struct CoolingTest {
 }; // dummy type to allow compile-type polymorphism via template specialization
 
-constexpr double m_H = quokka::hydrogen_mass_cgs;
+constexpr double m_H = C::m_u;
 constexpr double seconds_in_year = 3.154e7;
 
 template <> struct quokka::EOS_Traits<CoolingTest> {
 	static constexpr double gamma = 5. / 3.; // default value
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<CoolingTest> {
@@ -131,7 +131,7 @@ template <> void RadhydroSimulation<CoolingTest>::setInitialConditionsOnGrid(quo
 		Real xmom = 0;
 		Real ymom = 0;
 		Real zmom = 0;
-		Real const P = 4.0e4 * quokka::boltzmann_constant_cgs; // erg cm^-3
+		Real const P = 4.0e4 * C::k_B; // erg cm^-3
 		Real Eint = (quokka::EOS_Traits<CoolingTest>::gamma - 1.) * P;
 
 		Real const Egas = RadSystem<CoolingTest>::ComputeEgasFromEint(rho, xmom, ymom, zmom, Eint);
diff --git a/src/EOS.hpp b/src/EOS.hpp
index d324af982..7dbf8f75d 100644
--- a/src/EOS.hpp
+++ b/src/EOS.hpp
@@ -14,6 +14,7 @@
 #include "AMReX_GpuQualifiers.H"
 #include "AMReX_REAL.H"
 #include "physics_info.hpp"
+#include <AMReX_Print.H>
 
 #include "burn_type.H"
 #include "eos.H"
@@ -25,8 +26,6 @@
 
 namespace quokka
 {
-static constexpr double boltzmann_constant_cgs = 1.380658e-16; // cgs
-static constexpr double hydrogen_mass_cgs = 1.6726231e-24;     // cgs
 
 // specify default values for ideal gamma-law EOS
 //
@@ -34,7 +33,7 @@ template <typename problem_t> struct EOS_Traits {
 	static constexpr double gamma = 5. / 3.;     // default value
 	static constexpr double cs_isothermal = NAN; // only used when gamma = 1
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <typename problem_t> class EOS
@@ -84,8 +83,13 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeTgasFromEin
 	Tgas = chemstate.T;
 #else
 	if constexpr (gamma_ != 1.0) {
-		const amrex::Real c_v = boltzmann_constant_ / (mean_molecular_weight_ * (gamma_ - 1.0));
-		Tgas = Eint / (rho * c_v);
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.e = Eint / rho;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_re, estate);
+		// scale returned temperature in case boltzmann constant is dimensionless
+		Tgas = estate.T * C::k_B / boltzmann_constant_;
 	}
 #endif
 	return Tgas;
@@ -121,8 +125,12 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeEintFromTga
 	Eint = chemstate.e * chemstate.rho;
 #else
 	if constexpr (gamma_ != 1.0) {
-		const amrex::Real c_v = boltzmann_constant_ / (mean_molecular_weight_ * (gamma_ - 1.0));
-		Eint = rho * c_v * Tgas;
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.T = Tgas;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_rt, estate);
+		Eint = estate.e * rho * boltzmann_constant_ / C::k_B;
 	}
 #endif
 	return Eint;
@@ -130,7 +138,7 @@ AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto EOS<problem_t>::ComputeEintFromTga
 
 template <typename problem_t>
 AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE auto
-EOS<problem_t>::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Real /*Tgas*/,
+EOS<problem_t>::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Real Tgas,
 					  const std::optional<amrex::GpuArray<amrex::Real, nmscalars_>> massScalars) -> amrex::Real
 {
 	// compute derivative of internal energy w/r/t temperature
@@ -157,8 +165,12 @@ EOS<problem_t>::ComputeEintTempDerivative(const amrex::Real rho, const amrex::Re
 	dEint_dT = chemstate.dedT * chemstate.rho;
 #else
 	if constexpr (gamma_ != 1.0) {
-		const amrex::Real c_v = boltzmann_constant_ / (mean_molecular_weight_ * (gamma_ - 1.0));
-		dEint_dT = rho * c_v;
+		chem_eos_t estate;
+		estate.rho = rho;
+		estate.T = Tgas;
+		estate.mu = mean_molecular_weight_ / C::m_u;
+		eos(eos_input_rt, estate);
+		dEint_dT = estate.dedT * rho * boltzmann_constant_ / C::k_B;
 	}
 #endif
 	return dEint_dT;
diff --git a/src/FCQuantities/test_fc_quantities.cpp b/src/FCQuantities/test_fc_quantities.cpp
index bed07a180..38291fe22 100644
--- a/src/FCQuantities/test_fc_quantities.cpp
+++ b/src/FCQuantities/test_fc_quantities.cpp
@@ -27,8 +27,8 @@ struct FCQuantities {
 
 template <> struct quokka::EOS_Traits<FCQuantities> {
 	static constexpr double gamma = 5. / 3.;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<FCQuantities> {
diff --git a/src/HydroBlast2D/test_hydro2d_blast.cpp b/src/HydroBlast2D/test_hydro2d_blast.cpp
index 8ed15ee54..e8a3988e4 100644
--- a/src/HydroBlast2D/test_hydro2d_blast.cpp
+++ b/src/HydroBlast2D/test_hydro2d_blast.cpp
@@ -28,7 +28,7 @@ struct BlastProblem {
 template <> struct quokka::EOS_Traits<BlastProblem> {
 	static constexpr double gamma = 5. / 3.;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<BlastProblem> {
diff --git a/src/HydroBlast3D/test_hydro3d_blast.cpp b/src/HydroBlast3D/test_hydro3d_blast.cpp
index 42222dada..9a5e698a6 100644
--- a/src/HydroBlast3D/test_hydro3d_blast.cpp
+++ b/src/HydroBlast3D/test_hydro3d_blast.cpp
@@ -35,7 +35,7 @@ bool test_passes = false; // if one of the energy checks fails, set to false
 template <> struct quokka::EOS_Traits<SedovProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<SedovProblem> {
diff --git a/src/HydroContact/test_hydro_contact.cpp b/src/HydroContact/test_hydro_contact.cpp
index 411e5afdd..5576b7af3 100644
--- a/src/HydroContact/test_hydro_contact.cpp
+++ b/src/HydroContact/test_hydro_contact.cpp
@@ -23,8 +23,8 @@ struct ContactProblem {
 
 template <> struct quokka::EOS_Traits<ContactProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ContactProblem> {
diff --git a/src/HydroHighMach/test_hydro_highmach.cpp b/src/HydroHighMach/test_hydro_highmach.cpp
index a621e9453..82075656a 100644
--- a/src/HydroHighMach/test_hydro_highmach.cpp
+++ b/src/HydroHighMach/test_hydro_highmach.cpp
@@ -31,7 +31,7 @@ struct HighMachProblem {
 template <> struct quokka::EOS_Traits<HighMachProblem> {
 	static constexpr double gamma = 5. / 3.;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<HighMachProblem> {
diff --git a/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp b/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp
index 1885767c2..f561fe6f4 100644
--- a/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp
+++ b/src/HydroKelvinHelmholz/test_hydro2d_kh.cpp
@@ -25,7 +25,7 @@ struct KelvinHelmholzProblem {
 template <> struct quokka::EOS_Traits<KelvinHelmholzProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<KelvinHelmholzProblem> {
diff --git a/src/HydroLeblanc/test_hydro_leblanc.cpp b/src/HydroLeblanc/test_hydro_leblanc.cpp
index 9ce017f9a..7a848bb98 100644
--- a/src/HydroLeblanc/test_hydro_leblanc.cpp
+++ b/src/HydroLeblanc/test_hydro_leblanc.cpp
@@ -30,8 +30,8 @@ struct ShocktubeProblem {
 
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = (5. / 3.);
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
diff --git a/src/HydroQuirk/test_quirk.cpp b/src/HydroQuirk/test_quirk.cpp
index bc3d4ec6f..7ba3ac4b0 100644
--- a/src/HydroQuirk/test_quirk.cpp
+++ b/src/HydroQuirk/test_quirk.cpp
@@ -44,7 +44,7 @@ struct QuirkProblem {
 template <> struct quokka::EOS_Traits<QuirkProblem> {
 	static constexpr double gamma = 5. / 3.;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<QuirkProblem> {
diff --git a/src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp b/src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp
index 7dbbe13f4..05b1e5ae6 100644
--- a/src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp
+++ b/src/HydroRichtmeyerMeshkov/test_hydro2d_rm.cpp
@@ -24,7 +24,7 @@ struct RichtmeyerMeshkovProblem {
 template <> struct quokka::EOS_Traits<RichtmeyerMeshkovProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<RichtmeyerMeshkovProblem> {
diff --git a/src/HydroSMS/test_hydro_sms.cpp b/src/HydroSMS/test_hydro_sms.cpp
index b6495394a..ecfe6e42d 100644
--- a/src/HydroSMS/test_hydro_sms.cpp
+++ b/src/HydroSMS/test_hydro_sms.cpp
@@ -23,8 +23,8 @@ struct ShocktubeProblem {
 
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
diff --git a/src/HydroShocktube/test_hydro_shocktube.cpp b/src/HydroShocktube/test_hydro_shocktube.cpp
index 6211d0612..a1e6036f2 100644
--- a/src/HydroShocktube/test_hydro_shocktube.cpp
+++ b/src/HydroShocktube/test_hydro_shocktube.cpp
@@ -29,7 +29,7 @@ struct ShocktubeProblem {
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
diff --git a/src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp b/src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp
index 3a39ad0c8..fa7f921f1 100644
--- a/src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp
+++ b/src/HydroShocktubeCMA/test_hydro_shocktube_cma.cpp
@@ -34,7 +34,7 @@ template <> struct SimulationData<ShocktubeProblem> {
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
diff --git a/src/HydroShuOsher/test_hydro_shuosher.cpp b/src/HydroShuOsher/test_hydro_shuosher.cpp
index 5201c8ec3..3ce5b26c5 100644
--- a/src/HydroShuOsher/test_hydro_shuosher.cpp
+++ b/src/HydroShuOsher/test_hydro_shuosher.cpp
@@ -25,8 +25,8 @@ struct ShocktubeProblem {
 
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
diff --git a/src/HydroVacuum/test_hydro_vacuum.cpp b/src/HydroVacuum/test_hydro_vacuum.cpp
index a30dffdb6..42b8fcf99 100644
--- a/src/HydroVacuum/test_hydro_vacuum.cpp
+++ b/src/HydroVacuum/test_hydro_vacuum.cpp
@@ -28,7 +28,7 @@ struct ShocktubeProblem {
 template <> struct quokka::EOS_Traits<ShocktubeProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShocktubeProblem> {
diff --git a/src/HydroWave/test_hydro_wave.cpp b/src/HydroWave/test_hydro_wave.cpp
index 654501929..a69cb1329 100644
--- a/src/HydroWave/test_hydro_wave.cpp
+++ b/src/HydroWave/test_hydro_wave.cpp
@@ -23,8 +23,8 @@ struct WaveProblem {
 
 template <> struct quokka::EOS_Traits<WaveProblem> {
 	static constexpr double gamma = 5. / 3.;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<WaveProblem> {
diff --git a/src/ODEIntegration/test_ode.cpp b/src/ODEIntegration/test_ode.cpp
index c674ec1af..e5841b4b2 100644
--- a/src/ODEIntegration/test_ode.cpp
+++ b/src/ODEIntegration/test_ode.cpp
@@ -8,16 +8,18 @@
 ///
 
 #include "test_ode.hpp"
+#include "eos.H"
+#include "extern_parameters.H"
 #include "radiation_system.hpp"
 
 using amrex::Real;
 
-constexpr double Tgas0 = 6000.;				  // K
-constexpr double rho0 = 0.01 * quokka::hydrogen_mass_cgs; // g cm^-3
+constexpr double Tgas0 = 6000.;			  // K
+constexpr double rho0 = 0.01 * (C::m_p + C::m_e); // g cm^-3
 
 template <> struct quokka::EOS_Traits<ODETest> {
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };
 
@@ -30,7 +32,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto cooling_function(Real const rho, R
 	// use fitting function from Koyama & Inutsuka (2002)
 	Real gamma_heat = 2.0e-26;
 	Real lambda_cool = gamma_heat * (1.0e7 * std::exp(-114800. / (T + 1000.)) + 14. * std::sqrt(T) * std::exp(-92. / T));
-	Real rho_over_mh = rho / quokka::hydrogen_mass_cgs;
+	Real rho_over_mh = rho / (C::m_p + C::m_e);
 	Real cooling_source_term = rho_over_mh * gamma_heat - (rho_over_mh * rho_over_mh) * lambda_cool;
 	return cooling_source_term;
 }
@@ -52,6 +54,12 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE auto user_rhs(Real /*t*/, quokka::valar
 
 auto problem_main() -> int
 {
+	// initialize EOS
+	init_extern_parameters();
+	Real small_temp = 1e-10;
+	Real small_dens = 1e-100;
+	eos_init(small_temp, small_dens);
+
 	// set up initial conditions
 	const Real Eint0 = quokka::EOS<ODETest>::ComputeEintFromTgas(rho0, Tgas0);
 	const Real Edot0 = cooling_function(rho0, Tgas0);
diff --git a/src/PassiveScalar/test_scalars.cpp b/src/PassiveScalar/test_scalars.cpp
index 2bea41404..7f52c9859 100644
--- a/src/PassiveScalar/test_scalars.cpp
+++ b/src/PassiveScalar/test_scalars.cpp
@@ -25,8 +25,8 @@ struct ScalarProblem {
 
 template <> struct quokka::EOS_Traits<ScalarProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ScalarProblem> {
diff --git a/src/PrimordialChem/test_primordial_chem.cpp b/src/PrimordialChem/test_primordial_chem.cpp
index 8d2d4b946..e944dd166 100644
--- a/src/PrimordialChem/test_primordial_chem.cpp
+++ b/src/PrimordialChem/test_primordial_chem.cpp
@@ -40,13 +40,6 @@ using amrex::Real;
 struct PrimordialChemTest {
 }; // dummy type to allow compile-type polymorphism via template specialization
 
-// Currently, microphysics uses its own EOS, and this one below is used by hydro. Need to only have one EOS at some point.
-template <> struct quokka::EOS_Traits<PrimordialChemTest> {
-	static constexpr double gamma = 5. / 3.; // default value
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
-};
-
 template <> struct Physics_Traits<PrimordialChemTest> {
 	// cell-centred
 	static constexpr bool is_hydro_enabled = true;
diff --git a/src/RadBeam/test_radiation_beam.cpp b/src/RadBeam/test_radiation_beam.cpp
index ad1ea6148..1f4a4adf2 100644
--- a/src/RadBeam/test_radiation_beam.cpp
+++ b/src/RadBeam/test_radiation_beam.cpp
@@ -33,8 +33,8 @@ constexpr double a_rad = 7.5646e-15; // erg cm^-3 K^-4
 constexpr double c = 2.99792458e10;  // cm s^-1
 
 template <> struct quokka::EOS_Traits<BeamProblem> {
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };
 
diff --git a/src/RadForce/test_radiation_force.cpp b/src/RadForce/test_radiation_force.cpp
index bfa8f4cdc..265e28146 100644
--- a/src/RadForce/test_radiation_force.cpp
+++ b/src/RadForce/test_radiation_force.cpp
@@ -28,11 +28,11 @@
 struct TubeProblem {
 };
 
-constexpr double kappa0 = 5.0;				// cm^2 g^-1
-constexpr double mu = 2.33 * quokka::hydrogen_mass_cgs; // g
-constexpr double gamma_gas = 1.0;			// isothermal gas EOS
-constexpr double a0 = 0.2e5;				// cm s^-1
-constexpr double tau = 1.0e-6;				// optical depth (dimensionless)
+constexpr double kappa0 = 5.0;	     // cm^2 g^-1
+constexpr double mu = 2.33 * C::m_u; // g
+constexpr double gamma_gas = 1.0;    // isothermal gas EOS
+constexpr double a0 = 0.2e5;	     // cm s^-1
+constexpr double tau = 1.0e-6;	     // optical depth (dimensionless)
 
 constexpr double rho0 = 1.0e5 * mu; // g cm^-3
 constexpr double Mach0 = 1.1;	    // Mach number at wind base
@@ -44,8 +44,8 @@ constexpr double g0 = kappa0 * Frad0 / c_light_cgs_;	 // cm s^{-2}
 constexpr double Lx = (a0 * a0) / g0;			 // cm
 
 template <> struct quokka::EOS_Traits<TubeProblem> {
-	static constexpr double mean_molecular_mass = mu;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = mu;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = gamma_gas;
 	static constexpr double cs_isothermal = a0; // only used when gamma = 1
 };
diff --git a/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp b/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp
index d805efca0..9e99670ae 100644
--- a/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp
+++ b/src/RadMarshakAsymptotic/test_radiation_marshak_asymptotic.cpp
@@ -23,11 +23,11 @@ constexpr double T_initial = 300.;	      // K
 
 // constexpr double kelvin_to_eV = 8.617385e-5;
 constexpr double a_rad = radiation_constant_cgs_;
-constexpr double c_v = (quokka::boltzmann_constant_cgs / quokka::hydrogen_mass_cgs) / (5. / 3. - 1.);
+constexpr double c_v = (C::k_B / C::m_u) / (5. / 3. - 1.);
 
 template <> struct quokka::EOS_Traits<SuOlsonProblemCgs> {
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };
 
diff --git a/src/RadMarshakCGS/test_radiation_marshak_cgs.cpp b/src/RadMarshakCGS/test_radiation_marshak_cgs.cpp
index 8309cc633..39cdbf0ad 100644
--- a/src/RadMarshakCGS/test_radiation_marshak_cgs.cpp
+++ b/src/RadMarshakCGS/test_radiation_marshak_cgs.cpp
@@ -33,8 +33,8 @@ constexpr double alpha_SuOlson = 4.0 * a_rad / eps_SuOlson;
 constexpr double T_initial = 1.0e4; // K
 
 template <> struct quokka::EOS_Traits<SuOlsonProblemCgs> {
-	static constexpr double mean_molecular_mass = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_mass = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };
 
diff --git a/src/RadMatterCoupling/test_radiation_matter_coupling.cpp b/src/RadMatterCoupling/test_radiation_matter_coupling.cpp
index 515ab9e4b..3e554b767 100644
--- a/src/RadMatterCoupling/test_radiation_matter_coupling.cpp
+++ b/src/RadMatterCoupling/test_radiation_matter_coupling.cpp
@@ -31,8 +31,8 @@ template <> struct SimulationData<CouplingProblem> {
 };
 
 template <> struct quokka::EOS_Traits<CouplingProblem> {
-	static constexpr double mean_molecular_mass = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_mass = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };
 
diff --git a/src/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp b/src/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp
index b99dc52b6..3c7f2e97c 100644
--- a/src/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp
+++ b/src/RadMatterCouplingRSLA/test_radiation_matter_coupling_rsla.cpp
@@ -33,8 +33,8 @@ template <> struct SimulationData<CouplingProblem> {
 };
 
 template <> struct quokka::EOS_Traits<CouplingProblem> {
-	static constexpr double mean_molecular_mass = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_mass = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };
 
diff --git a/src/RadPulse/test_radiation_pulse.cpp b/src/RadPulse/test_radiation_pulse.cpp
index 49bd9ce92..e51da1d32 100644
--- a/src/RadPulse/test_radiation_pulse.cpp
+++ b/src/RadPulse/test_radiation_pulse.cpp
@@ -9,6 +9,7 @@
 
 #include "test_radiation_pulse.hpp"
 #include "AMReX_BC_TYPES.H"
+#include "AMReX_Print.H"
 #include "RadhydroSimulation.hpp"
 #include "fextract.hpp"
 
diff --git a/src/RadShadow/test_radiation_shadow.cpp b/src/RadShadow/test_radiation_shadow.cpp
index d414a0b29..1d11cb5e9 100644
--- a/src/RadShadow/test_radiation_shadow.cpp
+++ b/src/RadShadow/test_radiation_shadow.cpp
@@ -34,8 +34,8 @@ constexpr double a_rad = 7.5646e-15; // erg cm^-3 K^-4
 constexpr double c = 2.99792458e10;  // cm s^-1
 
 template <> struct quokka::EOS_Traits<ShadowProblem> {
-	static constexpr double mean_molecular_weight = 10. * quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = 10. * C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };
 
diff --git a/src/RadTophat/test_radiation_tophat.cpp b/src/RadTophat/test_radiation_tophat.cpp
index 613289074..5aae4dae4 100644
--- a/src/RadTophat/test_radiation_tophat.cpp
+++ b/src/RadTophat/test_radiation_tophat.cpp
@@ -38,8 +38,8 @@ constexpr double a_rad = 7.5646e-15; // erg cm^-3 K^-4
 constexpr double c = 2.99792458e10;  // cm s^-1
 
 template <> struct quokka::EOS_Traits<TophatProblem> {
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = 5. / 3.;
 };
 
diff --git a/src/RadTube/test_radiation_tube.cpp b/src/RadTube/test_radiation_tube.cpp
index c90f74bd4..13ddae128 100644
--- a/src/RadTube/test_radiation_tube.cpp
+++ b/src/RadTube/test_radiation_tube.cpp
@@ -25,8 +25,8 @@
 struct TubeProblem {
 };
 
-constexpr double kappa0 = 100.;				// cm^2 g^-1
-constexpr double mu = 2.33 * quokka::hydrogen_mass_cgs; // g
+constexpr double kappa0 = 100.;	     // cm^2 g^-1
+constexpr double mu = 2.33 * C::m_u; // g
 constexpr double gamma_gas = 5. / 3.;
 
 constexpr double rho0 = 1.0;		    // g cm^-3
@@ -38,7 +38,7 @@ constexpr double a0 = 4.0295519855200705e7; // cm s^-1
 
 template <> struct quokka::EOS_Traits<TubeProblem> {
 	static constexpr double mean_molecular_weight = mu;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 	static constexpr double gamma = gamma_gas;
 };
 
@@ -183,7 +183,7 @@ AMRSimulation<TubeProblem>::setCustomBoundaryConditions(const amrex::IntVect &iv
 		consVar(i, j, k, RadSystem<TubeProblem>::x2RadFlux_index) = 0.;
 		consVar(i, j, k, RadSystem<TubeProblem>::x3RadFlux_index) = 0.;
 
-		const double Egas = (quokka::boltzmann_constant_cgs / mu) * rho0 * T0 / (gamma_gas - 1.0);
+		const double Egas = (C::k_B / mu) * rho0 * T0 / (gamma_gas - 1.0);
 		const double x1Mom = consVar(lo[0], j, k, RadSystem<TubeProblem>::x1GasMomentum_index);
 		const double Ekin = 0.5 * (x1Mom * x1Mom) / rho0;
 		consVar(i, j, k, RadSystem<TubeProblem>::gasEnergy_index) = Egas + Ekin;
@@ -202,7 +202,7 @@ AMRSimulation<TubeProblem>::setCustomBoundaryConditions(const amrex::IntVect &iv
 		consVar(i, j, k, RadSystem<TubeProblem>::x2RadFlux_index) = 0;
 		consVar(i, j, k, RadSystem<TubeProblem>::x3RadFlux_index) = 0;
 
-		const double Egas = (quokka::boltzmann_constant_cgs / mu) * rho1 * T1 / (gamma_gas - 1.0);
+		const double Egas = (C::k_B / mu) * rho1 * T1 / (gamma_gas - 1.0);
 		const double x1Mom = consVar(hi[0], j, k, RadSystem<TubeProblem>::x1GasMomentum_index);
 		const double Ekin = 0.5 * (x1Mom * x1Mom) / rho1;
 		consVar(i, j, k, RadSystem<TubeProblem>::gasEnergy_index) = Egas + Ekin;
diff --git a/src/RadhydroShell/test_radhydro_shell.cpp b/src/RadhydroShell/test_radhydro_shell.cpp
index f6b2262fb..80e8487ca 100644
--- a/src/RadhydroShell/test_radhydro_shell.cpp
+++ b/src/RadhydroShell/test_radhydro_shell.cpp
@@ -36,12 +36,12 @@ struct ShellProblem {
 // if false, use octant symmetry
 constexpr bool simulate_full_box = true;
 
-constexpr double a_rad = 7.5646e-15;  // erg cm^-3 K^-4
-constexpr double c = 2.99792458e10;   // cm s^-1
-constexpr double a0 = 2.0e5;	      // ('reference' sound speed) [cm s^-1]
-constexpr double chat = 860. * a0;    // cm s^-1
-constexpr double k_B = 1.380658e-16;  // erg K^-1
-constexpr double m_H = 1.6726231e-24; // mass of hydrogen atom [g]
+constexpr double a_rad = 7.5646e-15; // erg cm^-3 K^-4
+constexpr double c = 2.99792458e10;  // cm s^-1
+constexpr double a0 = 2.0e5;	     // ('reference' sound speed) [cm s^-1]
+constexpr double chat = 860. * a0;   // cm s^-1
+constexpr double k_B = C::k_B;	     // erg K^-1
+constexpr double m_H = C::m_u;	     // atomic mass unit
 constexpr double gamma_gas = 5. / 3.;
 
 template <> struct quokka::EOS_Traits<ShellProblem> {
diff --git a/src/RadhydroShockCGS/test_radhydro_shock_cgs.cpp b/src/RadhydroShockCGS/test_radhydro_shock_cgs.cpp
index c0249b895..cafce3c0e 100644
--- a/src/RadhydroShockCGS/test_radhydro_shock_cgs.cpp
+++ b/src/RadhydroShockCGS/test_radhydro_shock_cgs.cpp
@@ -22,10 +22,9 @@ struct ShockProblem {
 // parameters taken from Section 9.5 of Skinner et al. (2019)
 // [The Astrophysical Journal Supplement Series, 241:7 (27pp), 2019 March]
 
-constexpr double a_rad = 7.5646e-15;  // erg cm^-3 K^-4
-constexpr double c = 2.99792458e10;   // cm s^-1
-constexpr double k_B = 1.380658e-16;  // erg K^-1
-constexpr double m_H = 1.6726231e-24; // mass of hydrogen atom [g]
+constexpr double a_rad = 7.5646e-15; // erg cm^-3 K^-4
+constexpr double c = 2.99792458e10;  // cm s^-1
+constexpr double k_B = C::k_B;	     // erg K^-1
 
 // constexpr double P0 = 1.0e-4;	// equal to P_0 in dimensionless units
 // constexpr double sigma_a = 1.0e6;	// absorption cross section
@@ -34,8 +33,7 @@ constexpr double c_s0 = 1.73e7; // adiabatic sound speed [cm s^-1]
 
 constexpr double kappa = 577.0; // "opacity" == rho*kappa [cm^-1] (!!)
 constexpr double gamma_gas = (5. / 3.);
-constexpr double mu = m_H;			       // mean molecular weight [grams]
-constexpr double c_v = k_B / (mu * (gamma_gas - 1.0)); // specific heat [erg g^-1 K^-1]
+constexpr double c_v = k_B / ((C::m_p + C::m_e) * (gamma_gas - 1.0)); // specific heat [erg g^-1 K^-1]
 
 constexpr double T0 = 2.18e6; // K
 constexpr double rho0 = 5.69; // g cm^-3
@@ -65,7 +63,7 @@ template <> struct RadSystem_Traits<ShockProblem> {
 };
 
 template <> struct quokka::EOS_Traits<ShockProblem> {
-	static constexpr double mean_molecular_weight = m_H;
+	static constexpr double mean_molecular_weight = C::m_p + C::m_e;
 	static constexpr double boltzmann_constant = k_B;
 	static constexpr double gamma = gamma_gas;
 };
diff --git a/src/RadhydroSimulation.hpp b/src/RadhydroSimulation.hpp
index e24f34caf..0ded9fc67 100644
--- a/src/RadhydroSimulation.hpp
+++ b/src/RadhydroSimulation.hpp
@@ -52,6 +52,8 @@
 #include "radiation_system.hpp"
 #include "simulation.hpp"
 
+#include "eos.H"
+
 // Simulation class should be initialized only once per program (i.e., is a singleton)
 template <typename problem_t> class RadhydroSimulation : public AMRSimulation<problem_t>
 {
@@ -127,6 +129,12 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 		defineComponentNames();
 		// read in runtime parameters
 		readParmParse();
+		// initialize Microphysics params
+		init_extern_parameters();
+		// initialize Microphysics EOS
+		amrex::Real small_temp = 1e-10;
+		amrex::Real small_dens = 1e-100;
+		eos_init(small_temp, small_dens);
 	}
 
 	explicit RadhydroSimulation(amrex::Vector<amrex::BCRec> &BCs_cc) : AMRSimulation<problem_t>(BCs_cc)
@@ -134,6 +142,12 @@ template <typename problem_t> class RadhydroSimulation : public AMRSimulation<pr
 		defineComponentNames();
 		// read in runtime parameters
 		readParmParse();
+		// initialize Microphysics params
+		init_extern_parameters();
+		// initialize Microphysics EOS
+		amrex::Real small_temp = 1e-10;
+		amrex::Real small_dens = 1e-100;
+		eos_init(small_temp, small_dens);
 	}
 
 	[[nodiscard]] static auto getScalarVariableNames() -> std::vector<std::string>;
diff --git a/src/RayleighTaylor2D/test_hydro2d_rt.cpp b/src/RayleighTaylor2D/test_hydro2d_rt.cpp
index 03bced418..6255e9b8f 100644
--- a/src/RayleighTaylor2D/test_hydro2d_rt.cpp
+++ b/src/RayleighTaylor2D/test_hydro2d_rt.cpp
@@ -25,7 +25,7 @@ struct RTProblem {
 template <> struct quokka::EOS_Traits<RTProblem> {
 	static constexpr double gamma = 1.4;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<RTProblem> {
diff --git a/src/RayleighTaylor3D/test_hydro3d_rt.cpp b/src/RayleighTaylor3D/test_hydro3d_rt.cpp
index cf4e66e0a..f05b814a3 100644
--- a/src/RayleighTaylor3D/test_hydro3d_rt.cpp
+++ b/src/RayleighTaylor3D/test_hydro3d_rt.cpp
@@ -24,8 +24,8 @@ struct RTProblem {
 
 template <> struct quokka::EOS_Traits<RTProblem> {
 	static constexpr double gamma = 1.4;
-	static constexpr double mean_molecular_weight = quokka::hydrogen_mass_cgs;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double mean_molecular_weight = C::m_u;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<RTProblem> {
diff --git a/src/ShockCloud/cloud.cpp b/src/ShockCloud/cloud.cpp
index c04a08615..a60435dd9 100644
--- a/src/ShockCloud/cloud.cpp
+++ b/src/ShockCloud/cloud.cpp
@@ -35,13 +35,12 @@ using amrex::Real;
 struct ShockCloud {
 }; // dummy type to allow compile-type polymorphism via template specialization
 
-constexpr double m_H = quokka::hydrogen_mass_cgs;
 constexpr double seconds_in_year = 3.154e7;
 
 template <> struct quokka::EOS_Traits<ShockCloud> {
 	static constexpr double gamma = 5. / 3.; // default value
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct Physics_Traits<ShockCloud> {
@@ -65,9 +64,9 @@ constexpr Real nH1 = 1.0e-1;		 // cm^-3
 constexpr Real R_cloud = 5.0 * 3.086e18; // cm [5 pc]
 constexpr Real M0 = 2.0;		 // Mach number of shock
 
-constexpr Real P0 = nH0 * Tgas0 * quokka::boltzmann_constant_cgs; // erg cm^-3
-constexpr Real rho0 = nH0 * m_H;				  // g cm^-3
-constexpr Real rho1 = nH1 * m_H;
+constexpr Real P0 = nH0 * Tgas0 * C::k_B;      // erg cm^-3
+constexpr Real rho0 = nH0 * (C::m_p + C::m_e); // g cm^-3
+constexpr Real rho1 = nH1 * (C::m_p + C::m_e);
 
 // perturbation parameters
 const int kmin = 0;
diff --git a/src/ShockCloud/test_cloudy.cpp b/src/ShockCloud/test_cloudy.cpp
index 2ce6560bc..336f8f095 100644
--- a/src/ShockCloud/test_cloudy.cpp
+++ b/src/ShockCloud/test_cloudy.cpp
@@ -50,15 +50,14 @@ auto problem_main() -> int
 	const Real T = ComputeTgasFromEgas(rho, Eint, quokka::EOS_Traits<ShockCloud>::gamma, tables);
 
 	const Real rhoH = rho * cloudy_H_mass_fraction;
-	const Real nH = rhoH / quokka::hydrogen_mass_cgs;
+	const Real nH = rhoH / (C::m_p + C::m_e);
 	const Real log_nH = std::log10(nH);
 
-	const Real C = (quokka::EOS_Traits<ShockCloud>::gamma - 1.) * Eint / (quokka::boltzmann_constant_cgs * (rho / quokka::hydrogen_mass_cgs));
+	const Real C = (quokka::EOS_Traits<ShockCloud>::gamma - 1.) * Eint / (C::k_B * (rho / (C::m_p + C::m_e)));
 	const Real mu = interpolate2d(log_nH, std::log10(T), tables.log_nH, tables.log_Tgas, tables.meanMolWeight);
 	const Real relerr = std::abs((C * mu - T) / T);
 
-	const Real n_e =
-	    (rho / quokka::hydrogen_mass_cgs) * (1.0 - mu * (X + Y / 4. + Z / mean_metals_A)) / (mu - (electron_mass_cgs / quokka::hydrogen_mass_cgs));
+	const Real n_e = (rho / (C::m_p + C::m_e)) * (1.0 - mu * (X + Y / 4. + Z / mean_metals_A)) / (mu - (electron_mass_cgs / (C::m_p + C::m_e)));
 
 	printf("\nrho = %.17e, Eint = %.17e, mu = %f, Tgas = %e, relerr = %e\n", rho, Eint, mu, T, relerr);
 	printf("n_e = %e, n_e/n_H = %e\n", n_e, n_e / nH);
diff --git a/src/SphericalCollapse/spherical_collapse.cpp b/src/SphericalCollapse/spherical_collapse.cpp
index a5273f0f2..253f289b5 100644
--- a/src/SphericalCollapse/spherical_collapse.cpp
+++ b/src/SphericalCollapse/spherical_collapse.cpp
@@ -29,7 +29,7 @@ struct CollapseProblem {
 template <> struct quokka::EOS_Traits<CollapseProblem> {
 	static constexpr double gamma = 5. / 3.;
 	static constexpr double mean_molecular_weight = NAN;
-	static constexpr double boltzmann_constant = quokka::boltzmann_constant_cgs;
+	static constexpr double boltzmann_constant = C::k_B;
 };
 
 template <> struct HydroSystem_Traits<CollapseProblem> {