diff --git a/src/allvars.cxx b/src/allvars.cxx
index 5de1a9bd..c8f78a8e 100644
--- a/src/allvars.cxx
+++ b/src/allvars.cxx
@@ -1811,6 +1811,9 @@ void PropDataHeader::declare_all_datasets(const Options &opt)
 #ifdef BHON
     declare_dataset<unsigned long>("n_bh");
     declare_dataset("Mass_bh", MASS);
+    //save ID most bound BH and its position
+    declare_dataset<long>("ID_mbp_bh");
+    declare_XYZ_datasets("cmbp_bh", LENGTH);
 #endif // BHON
 
 #ifdef HIGHRES
diff --git a/src/allvars.h b/src/allvars.h
index 1ed6d87f..415bdf48 100644
--- a/src/allvars.h
+++ b/src/allvars.h
@@ -1621,6 +1621,10 @@ struct PropData
     ///mean accretion rate, metallicty
     Double_t acc_bh, acc_bh_mostmassive;
 
+    ///ID of most bound BH and its position
+    Int_t ibound_bh = 0;
+    Coordinate gposmbp_bh;
+
     ///\name blackhole aperture/radial profiles
     //@{
     vector<int> aperture_npart_bh;
@@ -2478,6 +2482,7 @@ struct PropData
 #endif
 #ifdef BHON
         M_bh*=opt.h;
+        gposmbp_bh*=opt.h/opt.a;
 #endif
 #ifdef HIGHRES
         M_interloper*=opt.h;
diff --git a/src/io.cxx b/src/io.cxx
index 3cb3d92a..02c47917 100644
--- a/src/io.cxx
+++ b/src/io.cxx
@@ -2177,6 +2177,15 @@ void WriteProperties(Options &opt, const Int_t ngroups, PropData *pdata){
 
         for (Int_t i=0;i<ngroups;i++) ((Double_t*)data)[i]=pdata[i+1].M_bh;
         Fhdf.write_dataset(opt, head.headerdatainfo[itemp],ng,data,head.hdfpredtypeinfo[itemp]);itemp++;
+
+        for (Int_t i=0;i<ngroups;i++) ((long long*)data)[i]=pdata[i+1].ibound_bh;
+        Fhdf.write_dataset(opt, head.headerdatainfo[itemp],ng,data,head.hdfpredtypeinfo[itemp]);itemp++;
+
+        for (int k=0;k<3;k++){
+            for (Int_t i=0;i<ngroups;i++) ((Double_t*)data)[i]=pdata[i+1].gposmbp_bh[k];
+            Fhdf.write_dataset(opt, head.headerdatainfo[itemp],ng,data,head.hdfpredtypeinfo[itemp]);itemp++;
+        }
+
 #endif
 #ifdef HIGHRES
         for (Int_t i=0;i<ngroups;i++) ((unsigned long*)data)[i]=pdata[i+1].n_interloper;
diff --git a/src/substructureproperties.cxx b/src/substructureproperties.cxx
index 90489ec2..f5da189a 100644
--- a/src/substructureproperties.cxx
+++ b/src/substructureproperties.cxx
@@ -1136,6 +1136,14 @@ private(EncMassSF,EncMassNSF,Krot_sf,Krot_nsf,Ekin_sf,Ekin_nsf)
                 mval = opt.MassValue;
 #endif
                 pdata[i].M_bh+=mval;
+
+		//save particle information if the most bound id has not been assigned.
+		//note that this works because particles has been previously ordered by distance to defined centre.
+		if(pdata[i].ibound_bh == 0){
+			pdata[i].ibound_bh = Pval->GetPID();
+			for (k=0;k<3;k++) pdata[i].gposmbp_bh[k] = Pval->GetPosition(k) + pdata[i].gposmbp[k];
+		}
+
             }
         }
 #endif
@@ -2073,6 +2081,14 @@ private(j,Pval,x,y,z,vx,vy,vz,jval,jzval,zdist,Rdist)
                 mval = opt.MassValue;
 #endif
                 pdata[i].M_bh+=mval;
+
+                //save particle information if the most bound id has not been assigned.
+                //note that this works because particles has been previously ordered by distance to defined centre.
+                if(pdata[i].ibound_bh == 0){
+                        pdata[i].ibound_bh = Pval->GetPID();
+                        for (k=0;k<3;k++) pdata[i].gposmbp_bh[k] = Pval->GetPosition(k) + pdata[i].gposmbp[k];
+                }
+
             }
         }
 #endif