diff --git a/src/plugins/output_genericio.cc b/src/plugins/output_genericio.cc
index b5ba88e..f77bef8 100644
--- a/src/plugins/output_genericio.cc
+++ b/src/plugins/output_genericio.cc
@@ -10,7 +10,7 @@ protected:
     real_t lunit_, vunit_, munit_;
     bool hacc_hydro_;
     float hacc_etamax_;
-    float hh_value_, rho_value_, mu_value_;
+    float hh_value_, rho_value_, mu_value_, uu_value_;
     std::vector<int64_t> ids;
     std::vector<float> xx, yy, zz;
     std::vector<float> vx, vy, vz;
@@ -20,22 +20,43 @@ protected:
 
 public:
     //! constructor
-    explicit genericio_output_plugin(config_file &cf, cosmology::calculator &cc)
-            : output_plugin(cf, cc, "GenericIO")
+    explicit genericio_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator> &pcc)
+            : output_plugin(cf, pcc, "GenericIO")
     {
-        real_t astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
+        const double astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
+        const double hubble_param = pcc->cosmo_param_["h"];
         lunit_ = cf_.get_value<double>("setup", "BoxLength");
         vunit_ = lunit_/astart;
         hacc_hydro_ = cf_.get_value_safe<bool>("output", "GenericIO_HACCHydro", false);
         // initial smoothing length is mean particle seperation
         hh_value_ = lunit_ / cf_.get_value<float>("setup", "GridRes");
-        // neutral value for molecular weight. FIXME: account for ionization?
-        const float primordial_x = 0.75;
-        mu_value_ = 4.0 / (1.0 + 3.0*primordial_x);
-        
+
         double rhoc = 27.7519737 * 1e10; // in h^2 M_sol / Mpc^3
-        rho_value_ = cc_.cosmo_param_["Omega_b"] * rhoc;
         munit_ = rhoc * std::pow(cf_.get_value<double>("setup", "BoxLength"), 3);
+
+        if(hacc_hydro_) {
+            const double omegab = pcc_->cosmo_param_["Omega_b"];
+            const double gamma  = cf_.get_value_safe<double>("cosmology", "gamma", 5.0 / 3.0);
+            const double YHe    = pcc_->cosmo_param_["YHe"];
+            const double Tcmb0  = pcc_->cosmo_param_["Tcmb"];
+
+            // compute gas internal energy
+            const double npol = (fabs(1.0 - gamma) > 1e-7) ? 1.0 / (gamma - 1.) : 1.0;
+            const double unitv = 1e5; // km/s to cm/s
+            const double adec = 1.0 / (160. * std::pow(omegab * hubble_param * hubble_param / 0.022, 2.0 / 5.0));
+            const double Tini = astart < adec ? Tcmb0 / astart : Tcmb0 / astart / astart * adec;
+            const double mu = (Tini > 1.e4) ? 4.0 / (8. - 5. * YHe) : 4.0 / (1. + 3. * (1. - YHe));
+            const double ceint = 1.3806e-16 / 1.6726e-24 * Tini * npol / mu / unitv / unitv;
+
+            uu_value_ = ceint/astart/astart;  // probably needs a scale factor correction (might be in physical units).
+            mu_value_ = mu;
+
+            music::ilog.Print("HACC : calculated redshift of decoupling: z_dec = %.2f", 1./adec - 1.);
+            music::ilog.Print("HACC : set initial gas temperature to %.2e K/mu", Tini / mu);
+            music::ilog.Print("HACC : set initial internal energy to %.2e km^2/s^2", ceint);
+
+            rho_value_ = omegab * rhoc;
+        }
     }
 
     output_type write_species_as(const cosmo_species &) const { return output_type::particles; }
@@ -65,17 +86,15 @@ public:
         vz.reserve(vz.size() + npart);
         ids.reserve(ids.size() + npart);
         
-        
-
         auto _pos = reinterpret_cast<const float*>(pc.get_pos32_ptr());
         auto _vel = reinterpret_cast<const float*>(pc.get_vel32_ptr());
         auto _ids = reinterpret_cast<const uint64_t*>(pc.get_ids64_ptr());
         auto _mass = reinterpret_cast<const float*>(pc.get_mass32_ptr());
         
         for(size_t i=0; i<npart; ++i) {
-            xx.push_back(_pos[3*i + 0]);
-            yy.push_back(_pos[3*i + 1]);
-            zz.push_back(_pos[3*i + 2]);
+            xx.push_back(fmod(_pos[3*i + 0]+lunit_, lunit_));
+            yy.push_back(fmod(_pos[3*i + 1]+lunit_, lunit_));
+            zz.push_back(fmod(_pos[3*i + 2]+lunit_, lunit_));
             vx.push_back(_vel[3*i + 0]);
             vy.push_back(_vel[3*i + 1]);
             vz.push_back(_vel[3*i + 2]);
@@ -100,13 +119,13 @@ public:
             std::fill(mask.begin() + prev_size, mask.end(), s == cosmo_species::baryon ? 1<<2 : 0);
 
             hh.resize(new_size);
-            std::fill(hh.begin() + prev_size, hh.end(), s == cosmo_species::baryon ? hh_value_ : 0.0f);
+            std::fill(hh.begin() + prev_size, hh.end(), hh_value_);
 
             uu.resize(new_size);
-            std::fill(uu.begin() + prev_size, uu.end(), s == cosmo_species::baryon ? 0.0f : 0.0f);
+            std::fill(uu.begin() + prev_size, uu.end(), s == cosmo_species::baryon ? uu_value_ : 0.0f);
 
             rho.resize(new_size);
-            std::fill(rho.begin() + prev_size, rho.end(), s == cosmo_species::baryon ? rho_value_ : 0.0f);
+            std::fill(rho.begin() + prev_size, rho.end(), rho_value_);
 
             mu.resize(new_size);
             std::fill(mu.begin() + prev_size, mu.end(), s == cosmo_species::baryon ? mu_value_ : 0.0f);
@@ -141,7 +160,7 @@ public:
 
 namespace
 {
-output_plugin_creator_concrete<genericio_output_plugin> creator("genericio");
+output_plugin_creator_concrete<genericio_output_plugin> creator1("genericio");
 } // namespace
 
 #endif // ENABLE_GENERICIO
\ No newline at end of file