SWIFTSIM · annadurrant · Aug 15, 2026 · Oct 9, 2025 · Oct 10, 2025 · Oct 16, 2025
diff --git a/.gitignore b/.gitignore
@@ -15,8 +15,10 @@ documentation/table.tex
 documentation/timestamp.tex
 documentation/variations_table.tex
 documentation/units.tex
+documentation/property_table.rst
 
 tests/FLAMINGO/test_parameters.yml
 tests/COLIBRE/test_parameters.yml
 tests/test_SO_radius_*.png
 test_data/*
+test_SO_radius*png
diff --git a/CITATION.cff b/CITATION.cff
@@ -0,0 +1,37 @@
+cff-version: 1.2.0
+message: "If you use this software, please cite it as below."
+authors:
+  - family-names: "McGibbon"
+    given-names: "Robert"
+  - family-names: "Helly"
+    given-names: "John"
+  - family-names: "Schaye"
+    given-names: "Joop"
+  - family-names: "Schaller"
+    given-names: "Matthieu"
+  - family-names: "Vandenbroucke"
+    given-names: "Bert"
+title: "SOAP: A Python Package for Calculating the Properties of Galaxies and Halos Formed in Cosmological Simulations"
+doi: "10.21105/joss.08252"
+date-released: "2025-07-01"
+preferred-citation:
+  type: article
+  authors:
+    - family-names: "McGibbon"
+      given-names: "Robert"
+    - family-names: "Helly"
+      given-names: "John"
+    - family-names: "Schaye"
+      given-names: "Joop"
+    - family-names: "Schaller"
+      given-names: "Matthieu"
+    - family-names: "Vandenbroucke"
+      given-names: "Bert"
+  doi: "10.21105/joss.08252"
+  journal: "The Journal of Open Source Software"
+  month: 7
+  year: 2025
+  title: "SOAP: A Python Package for Calculating the Properties of Galaxies and Halos Formed in Cosmological Simulations"
+  volume: 10
+  issue: 111
+  start: 8252
diff --git a/README.md b/README.md
@@ -1,13 +1,15 @@
 # SOAP: Spherical Overdensity and Aperture Processor
 
-[![DOI](https://joss.theoj.org/papers/10.21105/joss.08252/status.svg)](https://doi.org/10.21105/joss.08252)
+[![DOI](https://joss.theoj.org/papers/10.21105/joss.08252/status.svg)](https://doi.org/10.21105/joss.08252) 
 
 This repository contains programs which can be used to compute
 properties of halos in spherical apertures in [SWIFT](https://swift.strw.leidenuniv.nl/) snapshots.
 The resulting output halo catalogues can be read using the
 [swiftsimio](https://swiftsimio.readthedocs.io/en/latest/)
 python package.
 
+Please cite SOAP using the [JOSS paper](https://ui.adsabs.harvard.edu/abs/2025JOSS...10.8252M)
+
 ## Installation
 
 The code is written in python and uses mpi4py for parallelism.
@@ -19,7 +21,7 @@ steps are recommended for install
 ```
 pip install mpi4py
 export HDF5_MPI="ON"; export CC=mpicc; pip install --no-binary=h5py h5py
-pip install git+https://github.com/SWIFTSIM/SOAP.git@soap_runtime
+pip install git+https://github.com/SWIFTSIM/SOAP.git
 ```
 
 ### Installation on COSMA
@@ -100,14 +102,66 @@ The optional `--max-ranks-reading` flag determines how many MPI ranks per node
 read the snapshot. This can be used to avoid overloading the file system. The
 default value is 32.
 
+### Selecting which subhalos to process
+
+By default SOAP calculates properties for every subhalo in the input catalogue.
+There are a number of flags which can be used to process only a subset of them.
+
+The `--centrals-only` flag discards satellites, so that only central subhalos are
+processed.
+
+Individual subhalos can be selected with the `--halo-indices` flag.
+This specifies the index of the required subhalos in the halo catalogue, which is
+the quantity written to `InputHalos/HaloCatalogueIndex`.
+
+For larger numbers of subhalos the indices can be listed in a text file, which is
+passed with the `--halo-indices-file` flag.
+The file must contain one index per line. Blank lines, and anything following a
+`#`, are ignored. Duplicate indices are discarded.
+
+### Command line arguments
+
+The arguments listed here are passed on the command line, and cannot be set in
+the parameter file.
+
+`SOAP/group_membership.py`:
+
+| Argument | Description |
+| --- | --- |
+| `config_file` | Name of the yaml parameter file. Required |
+| `--sim-name` | Name of the simulation to process |
+| `--snap-nr` | Snapshot number to process |
+
+`SOAP/compute_halo_properties.py`:
+
+| Argument | Default | Description |
+| --- | --- | --- |
+| `config_file` | | Name of the yaml parameter file. Required |
+| `--sim-name` | | Name of the simulation to process |
+| `--snap-nr` | | Snapshot number to process |
+| `--chunks` | 1 | Number of chunks to split the volume into. Should be at least the number of compute nodes |
+| `--dmo` | off | Run in dark matter only mode, skipping any hydro-only properties |
+| `--centrals-only` | off | Only process central halos, discarding satellites. See [Selecting which subhalos to process](#selecting-which-subhalos-to-process) |
+| `--halo-indices` | | Only process the listed halo indices. See [Selecting which subhalos to process](#selecting-which-subhalos-to-process) |
+| `--halo-indices-file` | | Only process the halo indices listed in the given file. See [Selecting which subhalos to process](#selecting-which-subhalos-to-process) |
+| `--max-halos` | 0 (all) | Only process the first N halos in the catalogue. See [Debugging](#debugging) |
+| `--record-halo-timings` | off | Record the time taken to process each halo. See [Timing](#timing) |
+| `--record-property-timings` | off | Record the time taken to calculate each property. This doubles the size of the output catalogue. See [Timing](#timing) |
+| `--reference-snapshot` | | Number of a snapshot which contains all particle types. Used to determine the datasets and units of any particle types which are missing from the snapshot being processed, e.g. stars or black holes at high redshift |
+| `--snipshot` / `--snapshot` | auto | Force snipshot or snapshot mode. By default this is determined from the value of `SelectOutput` in the snapshot header |
+| `--profile` | 0 | Run with profiling. See [Profiling](#profiling) |
+| `--max-ranks-reading` | 32 | Number of MPI ranks per node which read snapshot data. Can be reduced to avoid overloading the file system |
+| `--output-parameters` | | Where to write the parameters used by this run, in yaml format |
+
 ### Parameter files
 
 To run either of the programs a parameters file must be passed. This
 contains information including the input and output directories,
 the halo finder to use, which halo definitions to use, and
 which properties to calculate for each halo definition. A description
-of all possible fields, and a number of example parameter files
-can be found in the `parameters_files` directory.
+of all possible fields can be found in
+[`parameter_files/README.md`](parameter_files/README.md), alongside a number
+of example parameter files.
 
 ### Compression
 
@@ -116,13 +170,15 @@ The first is lossless compression via GZIP, the second is lossy compression.
 For the group membership files we only apply lossless compression. However,
 each property in the final SOAP catalogue has a lossy compression filter
 associated with it, which are set in `SOAP/property_table.py`. The script 
-`compression/compress_soap_catalogue.py` will apply both lossy and
+`SOAP/compression/compress_soap_catalogue.py` will apply both lossy and
 lossless compression to SOAP catalogues.
 
 ### Documentation
 
 A pdf describing the SOAP output can be generated. First run `SOAP/property_table.py` passing the parameter file used to run SOAP (to get the properties and halo types to include) and a snapshot (to get the units), e.g. `python SOAP/property_table.py parameter_files/COLIBRE_THERMAL.yml /cosma8/data/dp004/colibre/Runs/L0100N0752/Thermal/snapshots/colibre_0127/colibre_0127.hdf5`. This will generate a table containing all the properties which are enabled in the parameter file. To create the pdf run `cd documentation; pdflatex SOAP.tex; pdflatex SOAP.tex`. If you wish to see all possible properties then first run `python SOAP/property_table.py`, and then generate the pdf.
 
+The `property_table.py` script also generates the file `documentation/property_table.rst`, which is a version of the property table that can be useful for sphnix websites.
+
 ### Slurm scripts for running on COSMA
 
 The files in the `scripts` directory are made for running on cosma.
@@ -190,8 +246,9 @@ mpirun.
 The `-Werror` flag is useful for making pdb stop on warnings. E.g. division by
 zero in the halo property calculations will be caught.
 
-It is also possible to select individual halos to process with the `--halo-indices`
-flag. This specifies the index of the required halos in the halo catalogue. E.g.
+If SOAP crashes while processing a particular halo it will try to report the
+index of that halo, which can then be re-run on its own with the
+`--halo-indices` flag, e.g.
 ```
 python3 -Werror -m pdb ./compute_halo_properties.py --halo-indices 1 2 3 ...
 ```

diff --git a/SOAP/catalogue_readers/read_hbtplus.py b/SOAP/catalogue_readers/read_hbtplus.py
@@ -33,15 +33,15 @@ def read_hbtplus_groupnr(basename, read_potential_energies=False, registry=None)
     if comm_rank == 0:
         if os.path.exists(hbt_filename(basename, 0)):
             with h5py.File(hbt_filename(basename, 0), "r") as infile:
-                nr_files = int(infile["NumberOfFiles"][...])
+                nr_files = infile["NumberOfFiles"][0]
             sorted_file = False
         elif os.path.exists(basename):
             with h5py.File(basename, "r") as infile:
                 assert "Particles" in infile
             nr_files = 1
             sorted_file = True
         else:
-            print(f"No HBT files found for basename {basename}")
+            print(f"No HBT files found for basename {basename}", flush=True)
             comm.Abort()
     else:
         nr_files = None
@@ -241,9 +241,9 @@ def read_hbtplus_catalogue(
         have_units = False
         with h5py.File(filename, "r") as infile:
             if "Units" in infile:
-                LengthInMpch = float(infile["Units/LengthInMpch"][...])
-                MassInMsunh = float(infile["Units/MassInMsunh"][...])
-                VelInKmS = float(infile["Units/VelInKmS"][...])
+                LengthInMpch = infile["Units/LengthInMpch"][0]
+                MassInMsunh = infile["Units/MassInMsunh"][0]
+                VelInKmS = infile["Units/VelInKmS"][0]
                 have_units = True
         # Otherwise, will have to read the Parameters.log file
         if not (have_units):
@@ -264,7 +264,7 @@ def read_hbtplus_catalogue(
         MassInMsunh = None
         VelInKmS = None
         sorted_file = None
-    (LengthInMpch, MassInMsunh, VelInKmS) = comm.bcast(
+    LengthInMpch, MassInMsunh, VelInKmS = comm.bcast(
         (LengthInMpch, MassInMsunh, VelInKmS)
     )
     sorted_file = comm.bcast(sorted_file)

diff --git a/SOAP/catalogue_readers/read_subfind.py b/SOAP/catalogue_readers/read_subfind.py
@@ -201,7 +201,7 @@ def read_gadget4_catalogue(comm, basename, a_unit, registry, boxsize):
     length_conversion = (gadget_length_unit / swift_pmpc).to(unyt.dimensionless)
 
     # Get position in comoving Mpc, assuming input position from Gadget is comoving
-    cofp = data["Subhalo/SubhaloPos"] * length_conversion * swift_cmpc
+    cofp = (data["Subhalo/SubhaloPos"] * length_conversion * swift_cmpc) % boxsize
 
     # Store central halo flag
     is_central = np.where(data["Subhalo/SubhaloRankInGr"] == 0, 1, 0)
@@ -227,7 +227,7 @@ def read_gadget4_catalogue(comm, basename, a_unit, registry, boxsize):
 
     # Store initial search radius
     search_radius = (
-        data["Subhalo/SubhaloHalfmassRad"] * length_conversion * swift_pmpc
+        (4 * data["Subhalo/SubhaloHalfmassRad"]) * length_conversion * swift_pmpc
     )  # different units from cofm, not a typo!
 
     local_halo = {

diff --git a/SOAP/catalogue_readers/read_subfind_eagle.py b/SOAP/catalogue_readers/read_subfind_eagle.py
@@ -117,7 +117,8 @@ def read_subfind_catalogue(comm, basename, a_unit, registry, boxsize):
     )
 
     # Store initial search radius
-    search_radius = (5 * data["Subhalo/VmaxRadius"] / h) * swift_cmpc
+    search_radius_cmpc = np.minimum((5 * data["Subhalo/VmaxRadius"] / h), 5)
+    search_radius = search_radius_cmpc * swift_cmpc
 
     local_halo = {
         "cofp": cofp,

diff --git a/compression/README.md → SOAP/compression/README.md b/compression/README.md → SOAP/compression/README.md
diff --git a/compression/compress_soap_catalogue.py → SOAP/compression/compress_soap_catalogue.py b/compression/compress_soap_catalogue.py → SOAP/compression/compress_soap_catalogue.py
@@ -206,7 +206,7 @@ def assign_datasets(nr_files, nr_ranks, comm_rank):
             datasets = h5copy.dsets.copy()
 
         except Exception as e:
-            print(f"Error: {e}")
+            print(f"Error: {e}", flush=True)
             comm.Abort(1)
     else:
         tmp_dir = None

diff --git a/compression/create_empty_SOAP_catalogue.py → ...ompression/create_empty_SOAP_catalogue.py b/compression/create_empty_SOAP_catalogue.py → ...ompression/create_empty_SOAP_catalogue.py
@@ -115,7 +115,6 @@ def __call__(self, name, h5obj):
             elif name == "Parameters":
                 for attr in self.ifile[name].attrs:
                     self.ofile[name].attrs[attr] = self.ifile[name].attrs[attr]
-                self.ofile[name].attrs["halo_indices"] = np.array([], dtype="int64")
                 self.ofile[name].attrs["snapshot_nr"] = self.snapnum
             else:
                 for attr in self.ifile[name].attrs:

diff --git a/compression/extract_filters.py → SOAP/compression/extract_filters.py b/compression/extract_filters.py → SOAP/compression/extract_filters.py
diff --git a/compression/filters.yml → SOAP/compression/filters.yml b/compression/filters.yml → SOAP/compression/filters.yml