Compute kinetic and potential energy for DMO runs

SOAP/compute_halo_properties.py

@@ -294,6 +294,10 @@ def compute_halo_properties():
             )
 
     aperture_variations = parameter_file.get_halo_type_variations("ApertureProperties")
+    if args.dmo:
+        aperture_variations = parameter_file.remove_non_dmo_apertures(
+            "ApertureProperties", aperture_variations
+        )
 
     # Sort the aperture variations based on their radii, and create a list
     # of all apertures. This is required since we can skip some of the larger
@@ -405,6 +409,10 @@ def compute_halo_properties():
     projected_aperture_variations = parameter_file.get_halo_type_variations(
         "ProjectedApertureProperties"
     )
+    if args.dmo:
+        projected_aperture_variations = parameter_file.remove_non_dmo_apertures(
+            "ProjectedApertureProperties", projected_aperture_variations
+        )
     # Sort the aperture variations based on their radii, and create a list
     # of all apertures. This is required since we can skip some of the larger
     # apertures if all the particles were already included in the previous aperture
@@ -524,6 +532,7 @@ def compute_halo_properties():
         parameter_file.print_unregistered_properties(halo_prop_list, dmo=args.dmo)
         parameter_file.print_skipped_properties(halo_prop_list, dmo=args.dmo)
         parameter_file.print_optin_skipped_properties(halo_prop_list, dmo=args.dmo)
+        parameter_file.print_dmo_skipped_apertures()
         parameter_file.print_invalid_properties(halo_prop_list)
         parameter_file.print_variation_warnings()
         if not parameter_file.renclose_enabled():

SOAP/core/parameter_file.py

@@ -172,6 +172,10 @@ def __init__(
         # be calculated because the input files lack the datasets they need.
         self.uncomputable_properties = {}
 
+        # Apertures which are not calculated in a DMO run because their radius
+        # is defined by a non-DMO property. Stored as {variation: property}
+        self.dmo_skipped_apertures = {}
+
         # Names of the properties used by the filters defined in the parameter
         # file, generated on demand by _filter_property_names()
         self.filter_property_names = None
@@ -487,6 +491,47 @@ def print_optin_skipped_properties(
             for property in sorted(skipped):
                 print(f"  {property.ljust(40)}{skipped[property]}")
 
+    def remove_non_dmo_apertures(self, base_halo_type: str, variations: Dict) -> Dict:
+        """
+        Remove the variations whose aperture radius is defined by a property
+        which is not calculated in a DMO run. The removed variations are
+        recorded, and can be printed with print_dmo_skipped_apertures().
+
+        Parameters:
+         - base_halo_type: str
+           Halo type identifier in the parameter file, can be one of
+           ApertureProperties or ProjectedApertureProperties.
+         - variations: Dict
+           Dictionary of variations, as returned by get_halo_type_variations.
+
+        Returns the dictionary of variations which can be calculated.
+        """
+        kept = {}
+        for name, variation in variations.items():
+            property_name = variation.get("property")
+            if property_name is not None:
+                prop = _property_by_name(property_name.split("/")[-1])
+                if (prop is not None) and (not prop.dmo_property):
+                    self.dmo_skipped_apertures[f"{base_halo_type}/{name}"] = (
+                        property_name
+                    )
+                    continue
+            kept[name] = variation
+        return kept
+
+    def print_dmo_skipped_apertures(self) -> None:
+        """
+        Print a list of the apertures which are not calculated in a DMO run
+        because their radius is defined by a non-DMO property.
+        """
+        if not len(self.dmo_skipped_apertures):
+            return
+        print(
+            "Not computing the following apertures as they are defined by non-DMO properties:"
+        )
+        for name in sorted(self.dmo_skipped_apertures):
+            print(f"  {name.ljust(60)}{self.dmo_skipped_apertures[name]}")
+
     def print_uncomputable_properties(self) -> None:
         """
         Print a list of the properties which are enabled in the parameter file,

SOAP/property_table.py

@@ -1294,7 +1294,7 @@ class PropertyTable:
             unit="snap_mass*snap_length**2/snap_time**2",
             description="Total kinetic energy of the particles, relative to the centre of mass velocity.",
             lossy_compression_filter="FMantissa9",
-            dmo_property=False,
+            dmo_property=True,
             particle_properties=[
                 "PartType0/Masses",
                 "PartType0/Velocities",
@@ -1319,7 +1319,7 @@ class PropertyTable:
             unit="snap_mass*snap_length**2/snap_time**2",
             description="Total potential energy of the subhalo.",
             lossy_compression_filter="FMantissa9",
-            dmo_property=False,
+            dmo_property=True,
             particle_properties=[
                 "PartType0/SpecificPotentialEnergies",
                 "PartType0/Masses",