Update grackle

examples/GEAR/ZoomIn/zoom_in.yml

 # Define the system of units to use internally. 
 InternalUnitSystem:
-  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
-  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
+  UnitMass_in_cgs:     1.98841e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e24 # Mpc in centimeters
   UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
+Scheduler:
+  max_top_level_cells: 16
+  cell_sub_size_self_grav:   2048
+  cell_subdepth_diff_grav:   2
+  cell_split_size:           200
+  cell_extra_gparts:         100         # (Optional) Number of spare gparts per top-level allocated at rebuild time for on-the-fly creation.
+
+  
+# Parameters governing the time integration
+TimeIntegration:
+  dt_min:     1e-16 # The minimal time-step size of the simulation (in internal units).
+  dt_max:     0.1  # The maximal time-step size of the simulation (in internal units).
+  max_dt_RMS_factor: 0.25  # (Optional) Dimensionless factor for the maximal displacement allowed based on the RMS velocities.
+
 # Cosmological parameters
 Cosmology:
   h:              0.673        # Reduced Hubble constant
-  a_begin:        0.9873046739     # Initial scale-factor of the simulation
+  a_begin:        0.014084507042253521     # Initial scale-factor of the simulation
   a_end:          1.0           # Final scale factor of the simulation
   Omega_m:        0.315         # Matter density parameter
   Omega_lambda:   0.685         # Dark-energy density parameter
-  Omega_b:        0.0486        # Baryon density parameter
-  
-Scheduler:
-  max_top_level_cells:    8
-  
-# Parameters governing the time integration
-TimeIntegration:
-  time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1e-2  # The end time of the simulation (in internal units).
-  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
-  dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
-  
+  Omega_b:        0.0486     # Baryon density parameter
+  Omega_r:        0.            # (Optional) Radiation density parameter
+  w_0:            -1.0          # (Optional) Dark-energy equation-of-state parameter at z=0.
+  w_a:            0.            # (Optional) Dark-energy equation-of-state time evolution parameter.
+
 # Parameters governing the snapshots
 Snapshots:
-  basename:            snapshot # Common part of the name of output files
-  scale_factor_first:  0.987345  # Scale-factor of the first snaphot (cosmological run)
-  time_first:          0.01  # Time of the first output (non-cosmological run) (in internal units)
-  delta_time:          1.01  # Time difference between consecutive outputs (in internal units)
-  compression:         1
+  basename:            snap/h050 # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          1e-2  # Time difference between consecutive outputs (in internal units)
+  compression:         4
+  output_list_on:      1  # (Optional) Enable the output list
+  output_list:         snaplist.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
 
 # Parameters governing the conserved quantities statistics
 Statistics:
-  scale_factor_first:  0.987345 # Scale-factor of the first stat dump (cosmological run)
-  time_first:          0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
-  delta_time:          1.05 # Time between statistics output
+  scale_factor_first: 0.0141
+  delta_time:          1.02 # Time between statistics output
 
+# Parameters for the self-gravity scheme
 Gravity:
-  eta:                       0.002    # Constant dimensionless multiplier for time integration.
-  MAC:                       adaptive
-  theta_cr:                  0.7
-  epsilon_fmm:               0.001
-  use_tree_below_softening:  1
+  eta:                    0.05    # Constant dimensionless multiplier for time integration.
+  MAC:                           geometric  # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
+  epsilon_fmm:                   0.001     # Tolerance parameter for the adaptive multipole acceptance criterion.
+  theta_cr:                      0.7       # Opening angle for the purely gemoetric criterion.
   comoving_DM_softening:     0.1278 # Comoving softening length (in internal units).
   max_physical_DM_softening: 0.03365    # Physical softening length (in internal units).
   comoving_baryon_softening:     0.03365 # Comoving softening length (in internal units).
   max_physical_baryon_softening: 0.00673    # Physical softening length (in internal units).
   softening_ratio_background:    0.0285714      # Fraction of the mean inter-particle separation to use as Plummer-equivalent softening for the background DM particles.
   mesh_side_length:       128        # Number of cells along each axis for the periodic gravity mesh.
-
+  
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
   CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
-  minimal_temperature:   100      # (internal units)
+  minimal_temperature:   10.      # Kelvin
+  h_min_ratio:           0.0714285714285       # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified.
+  h_max:                 100     # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
 
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./h050.hdf5     # The file to read
-  periodic:   1
-  cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
-  cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
-
+  periodic:   1                     # Non-periodic BCs
+  cleanup_h_factors: 1              # Remove the h-factors inherited from Gadget
+  shift:    [0, 0, 0]   # Centre the box
+  cleanup_velocity_factors:    1    # (Optional) Clean up the scale-factors used in the definition of the velocity variable in the ICs (e.g. in Gadget files).
+  cleanup_smoothing_lengths:   1    # (Optional) Clean the values of the smoothing lengths that are read in to remove stupid values. Set to 1 to activate.
 
-# Cooling with Grackle 3.0
+# Cooling with Grackle 2.0
 GrackleCooling:
   cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
-  with_UV_background: 1                   # Enable or not the UV background
-  redshift: -1                           # Redshift to use (-1 means time based redshift)
-  with_metal_cooling: 1                   # Enable or not the metal cooling
-  provide_volumetric_heating_rates: 0      # (optional) User provide volumetric heating rates
-  provide_specific_heating_rates: 0        # (optional) User provide specific heating rates
-  self_shielding_method: 0                # (optional) Grackle (<= 3) or Gear self shielding method
-  max_steps: 10000                       # (optional) Max number of step when computing the initial composition
-  convergence_limit: 1e-2                # (optional) Convergence threshold (relative) for initial composition
-  thermal_time_myr:  5
-
-GEARChemistry:
-  initial_metallicity: 0.01295
+  with_UV_background: 1 # Enable or not the UV background
+  redshift: -1 # Redshift to use (-1 means time based redshift)
+  with_metal_cooling: 1 # Enable or not the metal cooling
+  provide_volumetric_heating_rates: 0 # User provide volumetric heating rates
+  provide_specific_heating_rates: 0 # User provide specific heating rates
+  self_shielding_method: -1 # Grackle (<= 3) or Gear self shielding method
+  self_shielding_threshold_atom_per_cm3: 0.007  # Required only with GEAR's self shielding. Density threshold of the self shielding
+  max_steps: 1000
+  convergence_limit: 1e-2
+  thermal_time_myr: 5
 
-GEARFeedback:
-  supernovae_energy_erg: 1e50
-  yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5
-  discrete_yields: 1
 
 GEARStarFormation:
   star_formation_efficiency: 0.01   # star formation efficiency (c_*)
@@ -94,4 +98,16 @@ GEARStarFormation:
   min_mass_frac: 0.5
 
 GEARPressureFloor:
-  jeans_factor: 10.       # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+  jeans_factor: 10
+
+GEARFeedback:
+  supernovae_energy_erg: 0.135e51
+  yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5
+  discrete_yields: 1
+
+GEARChemistry:
+  initial_metallicity: 0
+  scale_initial_metallicity: 0
+
+Restarts:
+  delta_hours:        1        # (Optional) decimal hours between dumps of restart files.

examples/parameter_example.yml

@@ -523,7 +523,7 @@ EAGLEFeedback:
 
 # GEAR feedback model
 GEARFeedback:
-  supernovae_energy_erg: 0.1e51                            # Energy released by a single supernovae.
+  supernovae_energy_erg: 0.135e51                            # Energy released by a single supernovae.
   yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5  # Table containing the yields.
   discrete_yields: 0                                       # Should we use discrete yields or the IMF integrated one?
 

src/chemistry/EAGLE/chemistry_io.h

@@ -50,11 +50,13 @@ INLINE static int chemistry_read_particles(struct part* parts,
  * @brief Specifies which particle fields to write to a dataset
  *
  * @param parts The particle array.
+ * @param xparts The extra particle array.
  * @param list The list of i/o properties to write.
  *
  * @return Returns the number of fields to write.
  */
 INLINE static int chemistry_write_particles(const struct part* parts,
+                                            const struct xpart* xparts,
                                             struct io_props* list) {
 
   /* List what we want to write */

src/chemistry/GEAR/chemistry.h

@@ -46,12 +46,17 @@
  * @param sp the new created star particle with its properties.
  */
 INLINE static void chemistry_copy_star_formation_properties(
-    const struct part* p, const struct xpart* xp, struct spart* sp) {
+    struct part* p, const struct xpart* xp, struct spart* sp) {
+
+  float mass = hydro_get_mass(p);
 
   /* Store the chemistry struct in the star particle */
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
     sp->chemistry_data.metal_mass_fraction[i] =
         p->chemistry_data.smoothed_metal_mass_fraction[i];
+
+    /* Remove the metals taken by the star. */
+    p->chemistry_data.metal_mass[i] *= mass / (mass + sp->mass);
   }
 }
 
@@ -174,6 +179,7 @@ __attribute__((always_inline)) INLINE static void chemistry_init_part(
     cpd->smoothed_metal_mass_fraction[i] = 0.f;
 
     /* Convert the total mass into mass fraction */
+    /* Now the metal mass is not available anymore */
     cpd->metal_mass_fraction[i] = cpd->metal_mass[i] / p->mass;
   }
 }
@@ -211,6 +217,7 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
     cpd->smoothed_metal_mass_fraction[i] *= factor;
 
     /* Convert the mass fraction into a total mass */
+    /* Now the metal mass fraction is not available anymore */
     cpd->metal_mass[i] = m * cpd->metal_mass_fraction[i];
   }
 }
@@ -370,7 +377,7 @@ __attribute__((always_inline)) INLINE static void chemistry_split_part(
  *
  * @param sp Pointer to the particle data.
  */
-__attribute__((always_inline)) INLINE static float
+__attribute__((always_inline)) INLINE static double
 chemistry_get_total_metal_mass_fraction_for_feedback(
     const struct spart* restrict sp) {
 
@@ -384,7 +391,7 @@ chemistry_get_total_metal_mass_fraction_for_feedback(
  *
  * @param sp Pointer to the particle data.
  */
-__attribute__((always_inline)) INLINE static float const*
+__attribute__((always_inline)) INLINE static double const*
 chemistry_get_metal_mass_fraction_for_feedback(
     const struct spart* restrict sp) {
 
@@ -397,7 +404,7 @@ chemistry_get_metal_mass_fraction_for_feedback(
  *
  * @param p Pointer to the particle data.
  */
-__attribute__((always_inline)) INLINE static float
+__attribute__((always_inline)) INLINE static double
 chemistry_get_total_metal_mass_fraction_for_cooling(
     const struct part* restrict p) {
 
@@ -411,7 +418,7 @@ chemistry_get_total_metal_mass_fraction_for_cooling(
  *
  * @param p Pointer to the particle data.
  */
-__attribute__((always_inline)) INLINE static float const*
+__attribute__((always_inline)) INLINE static double const*
 chemistry_get_metal_mass_fraction_for_cooling(const struct part* restrict p) {
 
   return p->chemistry_data.smoothed_metal_mass_fraction;
@@ -423,7 +430,7 @@ chemistry_get_metal_mass_fraction_for_cooling(const struct part* restrict p) {
  *
  * @param p Pointer to the particle data.
  */
-__attribute__((always_inline)) INLINE static float
+__attribute__((always_inline)) INLINE static double
 chemistry_get_total_metal_mass_fraction_for_star_formation(
     const struct part* restrict p) {
 
@@ -437,7 +444,7 @@ chemistry_get_total_metal_mass_fraction_for_star_formation(
  *
  * @param p Pointer to the particle data.
  */
-__attribute__((always_inline)) INLINE static float const*
+__attribute__((always_inline)) INLINE static double const*
 chemistry_get_metal_mass_fraction_for_star_formation(
     const struct part* restrict p) {
 

src/chemistry/GEAR/chemistry_io.h

@@ -42,33 +42,44 @@ INLINE static int chemistry_read_particles(struct part* parts,
 
   /* List what we want to read */
   list[0] = io_make_input_field(
-      "ElementAbundance", FLOAT, GEAR_CHEMISTRY_ELEMENT_COUNT, OPTIONAL,
+      "ElementAbundance", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT, OPTIONAL,
       UNIT_CONV_NO_UNITS, parts, chemistry_data.metal_mass_fraction);
 
   return 1;
 }
 
+INLINE static void convert_gas_metals(const struct engine* e,
+                                      const struct part* p,
+                                      const struct xpart* xp, double* ret) {
+
+  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
+    ret[0] = p->chemistry_data.metal_mass[i] / hydro_get_mass(p);
+  }
+}
+
 /**
  * @brief Specifies which particle fields to write to a dataset
  *
  * @param parts The particle array.
+ * @param xparts The extra particle array.
  * @param list The list of i/o properties to write.
  *
  * @return Returns the number of fields to write.
  */
 INLINE static int chemistry_write_particles(const struct part* parts,
+                                            const struct xpart* xparts,
                                             struct io_props* list) {
 
   /* List what we want to write */
   list[0] = io_make_output_field(
-      "SmoothedElementAbundances", FLOAT, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      "SmoothedElementAbundances", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
       UNIT_CONV_NO_UNITS, 0.f, parts,
       chemistry_data.smoothed_metal_mass_fraction,
       "Element abundances smoothed over the neighbors");
 
-  list[1] = io_make_output_field(
-      "ElementAbundances", FLOAT, GEAR_CHEMISTRY_ELEMENT_COUNT,
-      UNIT_CONV_NO_UNITS, 0.f, parts, chemistry_data.metal_mass_fraction,
+  list[1] = io_make_output_field_convert_part(
+      "ElementAbundances", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      UNIT_CONV_NO_UNITS, 0.f, parts, xparts, convert_gas_metals,
       "Mass fraction of each element");
 
   return 2;
@@ -87,7 +98,7 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
 
   /* List what we want to write */
   list[0] = io_make_output_field(
-      "ElementAbundances", FLOAT, GEAR_CHEMISTRY_ELEMENT_COUNT,
+      "ElementAbundances", DOUBLE, GEAR_CHEMISTRY_ELEMENT_COUNT,
       UNIT_CONV_NO_UNITS, 0.f, sparts, chemistry_data.metal_mass_fraction,
       "Mass fraction of each element");
 

src/chemistry/GEAR/chemistry_struct.h

@@ -25,7 +25,7 @@
 struct chemistry_global_data {
 
   /* Initial metallicity Z */
-  float initial_metallicities[GEAR_CHEMISTRY_ELEMENT_COUNT];
+  double initial_metallicities[GEAR_CHEMISTRY_ELEMENT_COUNT];
 };
 
 /**
@@ -33,14 +33,18 @@ struct chemistry_global_data {
  */
 struct chemistry_part_data {
 
-  /*! Fraction of the particle mass in a given element */
-  float metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
+  union {
+    /*! Fraction of the particle mass in a given element.
+      This field is available only during the density hydro loop. */
+    double metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
 
-  /*! Total mass of element in a particle */
-  float metal_mass[GEAR_CHEMISTRY_ELEMENT_COUNT];
+    /*! Total mass of element in a particle.
+      This field is available only outside the density hydro loop. */
+    double metal_mass[GEAR_CHEMISTRY_ELEMENT_COUNT];
+  };
 
   /*! Smoothed fraction of the particle mass in a given element */
-  float smoothed_metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
+  double smoothed_metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
 };
 
 /**
@@ -49,7 +53,7 @@ struct chemistry_part_data {
 struct chemistry_spart_data {
 
   /*! Fraction of the particle mass in a given element */
-  float metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
+  double metal_mass_fraction[GEAR_CHEMISTRY_ELEMENT_COUNT];
 };
 
 /**

src/chemistry/QLA/chemistry_io.h

@@ -42,11 +42,13 @@ INLINE static int chemistry_read_particles(struct part* parts,
  * @brief Specifies which particle fields to write to a dataset
  *
  * @param parts The particle array.
+ * @param xparts The extra particle array.
  * @param list The list of i/o properties to write.
  *
  * @return Returns the number of fields to write.
  */
 INLINE static int chemistry_write_particles(const struct part* parts,
+                                            const struct xpart* xparts,
                                             struct io_props* list) {
 
   /* update list according to hydro_io */

src/chemistry/none/chemistry_io.h

@@ -42,11 +42,13 @@ INLINE static int chemistry_read_particles(struct part* parts,
  * @brief Specifies which particle fields to write to a dataset
  *
  * @param parts The particle array.
+ * @param xparts The extra particle array.
  * @param list The list of i/o properties to write.
  *
  * @return Returns the number of fields to write.
  */
 INLINE static int chemistry_write_particles(const struct part* parts,
+                                            const struct xpart* xparts,
                                             struct io_props* list) {
 
   /* update list according to hydro_io */

src/common_io.c

@@ -2855,7 +2855,7 @@ int get_ptype_fields(const int ptype, struct io_props* list,
 
     case swift_type_gas:
       hydro_write_particles(NULL, NULL, list, &num_fields);
-      num_fields += chemistry_write_particles(NULL, list + num_fields);
+      num_fields += chemistry_write_particles(NULL, NULL, list + num_fields);
       num_fields +=
           cooling_write_particles(NULL, NULL, list + num_fields, NULL);
       num_fields += tracers_write_particles(NULL, NULL, list + num_fields,

src/cooling/grackle/cooling.c

@@ -284,6 +284,7 @@ void cooling_print_backend(const struct cooling_function_data* cooling) {
   if (cooling->self_shielding_method == -1) {
     message("Self Shelding density = %g", cooling->self_shielding_threshold);
   }
+  message("Thermal time = %g", cooling->thermal_time);
   message("Specific Heating Rates = %i",
           cooling->provide_specific_heating_rates);
   message("Volumetric Heating Rates = %i",

src/cooling/grackle/cooling_io.h

@@ -191,8 +191,8 @@ __attribute__((always_inline)) INLINE static void cooling_read_parameters(
       parameter_file, "GrackleCooling:convergence_limit", 1e-2);
 
   /* Thermal time */
-  cooling->thermal_time =
-      parser_get_param_float(parameter_file, "GrackleCooling:thermal_time_myr");
+  cooling->thermal_time = parser_get_param_double(
+      parameter_file, "GrackleCooling:thermal_time_myr");
   cooling->thermal_time *= phys_const->const_year * 1e6;
 }
 

src/cooling/grackle/cooling_struct.h

@@ -76,7 +76,7 @@ struct cooling_function_data {
   float omega;
 
   /*! Duration for switching off cooling after an event (e.g. supernovae) */
-  float thermal_time;
+  double thermal_time;
 };
 
 /**
@@ -93,7 +93,7 @@ struct cooling_xpart_data {
   float radiated_energy;
 
   /*! Last time the cooling was switch off */
-  float time_last_event;
+  double time_last_event;
 
 /* here all fractions are mass fraction */
 #if COOLING_GRACKLE_MODE >= 1

src/distributed_io.c

@@ -498,7 +498,8 @@ void write_output_distributed(struct engine* e,
           /* No inhibted particles: easy case */
           Nparticles = Ngas;
           hydro_write_particles(parts, xparts, list, &num_fields);
-          num_fields += chemistry_write_particles(parts, list + num_fields);
+          num_fields +=
+              chemistry_write_particles(parts, xparts, list + num_fields);
           if (with_cooling || with_temperature) {
             num_fields += cooling_write_particles(
                 parts, xparts, list + num_fields, e->cooling_func);
@@ -537,8 +538,8 @@ void write_output_distributed(struct engine* e,
           /* Select the fields to write */
           hydro_write_particles(parts_written, xparts_written, list,
                                 &num_fields);
-          num_fields +=
-              chemistry_write_particles(parts_written, list + num_fields);
+          num_fields += chemistry_write_particles(parts_written, xparts_written,
+                                                  list + num_fields);
           if (with_cooling || with_temperature) {
             num_fields +=
                 cooling_write_particles(parts_written, xparts_written,

src/feedback/GEAR/feedback_struct.h

@@ -57,7 +57,7 @@ struct feedback_spart_data {
   float mass_ejected;
 
   /*! Chemical composition of the mass ejected */
-  float metal_mass_ejected[GEAR_CHEMISTRY_ELEMENT_COUNT];
+  double metal_mass_ejected[GEAR_CHEMISTRY_ELEMENT_COUNT];
 };
 
 #endif /* SWIFT_FEEDBACK_STRUCT_GEAR_H */

src/line_of_sight.c

@@ -438,7 +438,7 @@ void write_los_hdf5_datasets(hid_t grp, const int j, const size_t N,
 
   /* Find all the gas output fields */
   hydro_write_particles(parts, xparts, list, &num_fields);
-  num_fields += chemistry_write_particles(parts, list + num_fields);
+  num_fields += chemistry_write_particles(parts, xparts, list + num_fields);
   if (with_cooling || with_temperature) {
     num_fields += cooling_write_particles(parts, xparts, list + num_fields,
                                           e->cooling_func);

src/parallel_io.c

@@ -1249,7 +1249,8 @@ void prepare_file(struct engine* e, const char* fileName,
 
       case swift_type_gas:
         hydro_write_particles(parts, xparts, list, &num_fields);
-        num_fields += chemistry_write_particles(parts, list + num_fields);
+        num_fields +=
+            chemistry_write_particles(parts, xparts, list + num_fields);
         if (with_cooling || with_temperature) {
           num_fields += cooling_write_particles(
               parts, xparts, list + num_fields, e->cooling_func);
@@ -1621,7 +1622,8 @@ void write_output_parallel(struct engine* e,
           /* No inhibted particles: easy case */
           Nparticles = Ngas;
           hydro_write_particles(parts, xparts, list, &num_fields);
-          num_fields += chemistry_write_particles(parts, list + num_fields);
+          num_fields +=
+              chemistry_write_particles(parts, xparts, list + num_fields);
           if (with_cooling || with_temperature) {
             num_fields += cooling_write_particles(
                 parts, xparts, list + num_fields, e->cooling_func);
@@ -1660,8 +1662,8 @@ void write_output_parallel(struct engine* e,
           /* Select the fields to write */
           hydro_write_particles(parts_written, xparts_written, list,
                                 &num_fields);
-          num_fields +=
-              chemistry_write_particles(parts_written, list + num_fields);
+          num_fields += chemistry_write_particles(parts_written, xparts_written,
+                                                  list + num_fields);
           if (with_cooling || with_temperature) {
             num_fields +=
                 cooling_write_particles(parts_written, xparts_written,

src/serial_io.c

@@ -1205,7 +1205,8 @@ void write_output_serial(struct engine* e,
               /* No inhibted particles: easy case */
               Nparticles = Ngas;
               hydro_write_particles(parts, xparts, list, &num_fields);
-              num_fields += chemistry_write_particles(parts, list + num_fields);
+              num_fields +=
+                  chemistry_write_particles(parts, xparts, list + num_fields);
               if (with_cooling || with_temperature) {
                 num_fields += cooling_write_particles(
                     parts, xparts, list + num_fields, e->cooling_func);
@@ -1244,8 +1245,8 @@ void write_output_serial(struct engine* e,
               /* Select the fields to write */
               hydro_write_particles(parts_written, xparts_written, list,
                                     &num_fields);
-              num_fields +=
-                  chemistry_write_particles(parts_written, list + num_fields);
+              num_fields += chemistry_write_particles(
+                  parts_written, xparts_written, list + num_fields);
               if (with_cooling || with_temperature) {
                 num_fields +=
                     cooling_write_particles(parts_written, xparts_written,

src/single_io.c

@@ -997,7 +997,8 @@ void write_output_single(struct engine* e,
           /* No inhibted particles: easy case */
           N = Ngas;
           hydro_write_particles(parts, xparts, list, &num_fields);
-          num_fields += chemistry_write_particles(parts, list + num_fields);
+          num_fields +=
+              chemistry_write_particles(parts, xparts, list + num_fields);
           if (with_cooling || with_temperature) {
             num_fields += cooling_write_particles(
                 parts, xparts, list + num_fields, e->cooling_func);
@@ -1036,8 +1037,8 @@ void write_output_single(struct engine* e,
           /* Select the fields to write */
           hydro_write_particles(parts_written, xparts_written, list,
                                 &num_fields);
-          num_fields +=
-              chemistry_write_particles(parts_written, list + num_fields);
+          num_fields += chemistry_write_particles(parts_written, xparts_written,
+                                                  list + num_fields);
           if (with_cooling || with_temperature) {
             num_fields +=
                 cooling_write_particles(parts_written, xparts_written,