diff --git a/src/pgen/cluster.cpp b/src/pgen/cluster.cpp
index ca1916e4..20eaf903 100644
--- a/src/pgen/cluster.cpp
+++ b/src/pgen/cluster.cpp
@@ -64,131 +64,6 @@ using namespace parthenon::package::prelude;
 using utils::few_modes_ft::FewModesFT;
 using utils::few_modes_ft_log::FewModesFTLog;
 
-<<<<<<< HEAD
-template <class EOS>
-void ApplyClusterClips(MeshData<Real> *md, const parthenon::SimTime &tm,
-                       const Real beta_dt, const EOS eos) {
-
-  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
-
-  // Apply clips -- ceilings on temperature, velocity, alfven velocity, and
-  // density floor -- within a radius of the AGN
-  const auto &dfloor = hydro_pkg->Param<Real>("cluster_dfloor");
-  const auto &eceil = hydro_pkg->Param<Real>("cluster_eceil");
-  const auto &vceil = hydro_pkg->Param<Real>("cluster_vceil");
-  const auto &vAceil = hydro_pkg->Param<Real>("cluster_vAceil");
-  const auto &clip_r = hydro_pkg->Param<Real>("cluster_clip_r");
-
-  if (clip_r > 0 && (dfloor > 0 || eceil < std::numeric_limits<Real>::infinity() ||
-                     vceil < std::numeric_limits<Real>::infinity() ||
-                     vAceil < std::numeric_limits<Real>::infinity())) {
-    // Grab some necessary variables
-    const auto &prim_pack = md->PackVariables(std::vector<std::string>{"prim"});
-    const auto &cons_pack = md->PackVariables(std::vector<std::string>{"cons"});
-    IndexRange ib = md->GetBlockData(0)->GetBoundsI(IndexDomain::interior);
-    IndexRange jb = md->GetBlockData(0)->GetBoundsJ(IndexDomain::interior);
-    IndexRange kb = md->GetBlockData(0)->GetBoundsK(IndexDomain::interior);
-    const auto nhydro = hydro_pkg->Param<int>("nhydro");
-    const auto nscalars = hydro_pkg->Param<int>("nscalars");
-
-    const Real clip_r2 = SQR(clip_r);
-    const Real vceil2 = SQR(vceil);
-    const Real vAceil2 = SQR(vAceil);
-    const Real gm1 = (hydro_pkg->Param<Real>("AdiabaticIndex") - 1.0);
-
-    parthenon::par_for(
-        DEFAULT_LOOP_PATTERN, "Cluster::ApplyClusterClips", parthenon::DevExecSpace(), 0,
-        cons_pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
-        KOKKOS_LAMBDA(const int &b, const int &k, const int &j, const int &i) {
-          auto &cons = cons_pack(b);
-          auto &prim = prim_pack(b);
-          const auto &coords = cons_pack.GetCoords(b);
-
-          const Real r2 =
-              SQR(coords.Xc<1>(i)) + SQR(coords.Xc<2>(j)) + SQR(coords.Xc<3>(k));
-
-          if (r2 < clip_r2) {
-            // Cell falls within clipping radius
-              
-            const int gks = 0;
-            const int gjs = 0;
-            const int gis = 0;
-              
-            eos.ConsToPrim(cons, prim, nhydro, nscalars, k, j, i); // Three last parameters are just passive here
-
-            if (dfloor > 0) {
-              const Real rho = prim(IDN, k, j, i);
-              if (rho < dfloor) {
-                cons(IDN, k, j, i) = dfloor;
-                prim(IDN, k, j, i) = dfloor;
-              }
-            }
-
-            if (vceil < std::numeric_limits<Real>::infinity()) {
-              // Apply velocity ceiling
-              const Real v2 = SQR(prim(IV1, k, j, i)) + SQR(prim(IV2, k, j, i)) +
-                              SQR(prim(IV3, k, j, i));
-              if (v2 > vceil2) {
-                // Fix the velocity to the velocity ceiling
-                const Real v = sqrt(v2);
-                cons(IM1, k, j, i) *= vceil / v;
-                cons(IM2, k, j, i) *= vceil / v;
-                cons(IM3, k, j, i) *= vceil / v;
-                prim(IV1, k, j, i) *= vceil / v;
-                prim(IV2, k, j, i) *= vceil / v;
-                prim(IV3, k, j, i) *= vceil / v;
-
-                // Remove kinetic energy
-                cons(IEN, k, j, i) -= 0.5 * prim(IDN, k, j, i) * (v2 - vceil2);
-              }
-            }
-
-            if (vAceil2 < std::numeric_limits<Real>::infinity()) {
-              // Apply Alfven velocity ceiling by raising density
-              const Real rho = prim(IDN, k, j, i);
-              const Real B2 = (SQR(prim(IB1, k, j, i)) + SQR(prim(IB2, k, j, i)) +
-                               SQR(prim(IB3, k, j, i)));
-
-              // compute Alfven mach number
-              const Real va2 = (B2 / rho);
-
-              if (va2 > vAceil2) {
-                // Increase the density to match the alfven velocity ceiling
-                const Real rho_new = std::sqrt(B2 / vAceil2);
-                cons(IDN, k, j, i) = rho_new;
-                prim(IDN, k, j, i) = rho_new;
-              }
-            }
-
-            if (eceil < std::numeric_limits<Real>::infinity()) {
-              // Apply  internal energy ceiling as a pressure ceiling
-              const Real internal_e = prim(IPR, k, j, i) / (gm1 * prim(IDN, k, j, i));
-              if (internal_e > eceil) {
-                cons(IEN, k, j, i) -= prim(IDN, k, j, i) * (internal_e - eceil);
-                prim(IPR, k, j, i) = gm1 * prim(IDN, k, j, i) * eceil;
-              }
-            }
-          }
-        });
-  }
-}
-
-void ApplyClusterClips(MeshData<Real> *md, const parthenon::SimTime &tm,
-                       const Real beta_dt) {
-  auto hydro_pkg = md->GetBlockData(0)->GetBlockPointer()->packages.Get("Hydro");
-  auto fluid = hydro_pkg->Param<Fluid>("fluid");
-  if (fluid == Fluid::euler) {
-    ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param<AdiabaticHydroEOS>("eos"));
-  } else if (fluid == Fluid::glmmhd) {
-    ApplyClusterClips(md, tm, beta_dt, hydro_pkg->Param<AdiabaticGLMMHDEOS>("eos"));
-  } else {
-    PARTHENON_FAIL("Cluster::ApplyClusterClips: Unknown EOS");
-  }
-}
-
-//void ClusterSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
-//                    const Real beta_dt) {
-
 void ClusterUnsplitSrcTerm(MeshData<Real> *md, const parthenon::SimTime &tm,
                            const Real beta_dt) {