diff --git a/src/include/const_accessor_simple.hxx b/src/include/const_accessor_simple.hxx
index e0c9b82c19..9e751f9876 100644
--- a/src/include/const_accessor_simple.hxx
+++ b/src/include/const_accessor_simple.hxx
@@ -9,7 +9,6 @@ struct ParticleProxySimple
 {
   using Particle = typename Mparticles::Particle;
   using real_t = typename Mparticles::real_t;
-  using Patch = typename Mparticles::Patch;
   using Real3 = Vec3<real_t>;
 
   ParticleProxySimple(Particle& prt, const Mparticles& mprts)
@@ -212,11 +211,8 @@ struct ConstAccessorSimple
 
   Patch operator[](int p) const { return {*this, p}; }
   const Mparticles& mprts() const { return mprts_; }
-  uint size(int p) const { return mprts_[p].size(); }
-  typename Mparticles::Patch::iterator data(int p) const
-  {
-    return mprts_[p].begin();
-  }
+  uint size(int p) const { return mprts_.size(p); }
+  typename Mparticles::iterator data(int p) const { return mprts_.begin(p); }
   const Grid_t& grid() const { return mprts_.grid(); }
 
 private:
@@ -237,8 +233,8 @@ struct AccessorSimple
   Patch operator[](int p) { return {*this, p}; }
   const Mparticles& mprts() const { return mprts_; }
   Mparticles& mprts() { return mprts_; }
-  uint size(int p) const { return mprts_[p].size(); }
-  typename Mparticles::Patch::iterator data(int p) { return mprts_[p].begin(); }
+  uint size(int p) const { return mprts_.size(p); }
+  typename Mparticles::iterator data(int p) { return mprts_.begin(p); }
   const Grid_t& grid() const { return mprts_.grid(); }
 
 private:
diff --git a/src/include/injector_simple.hxx b/src/include/injector_simple.hxx
index 0b5048ef08..ac4b49d034 100644
--- a/src/include/injector_simple.hxx
+++ b/src/include/injector_simple.hxx
@@ -32,7 +32,7 @@ struct InjectorSimple
         new_prt.kind,
         mprts_.uid_gen(),
         new_prt.tag};
-      mprts_[p_].push_back(prt);
+      mprts_.push_back(p_, prt);
     }
 
     void reweight(const psc::particle::Inject& new_prt)
diff --git a/src/include/particle_simple.hxx b/src/include/particle_simple.hxx
index 4f1b3e2107..6a504fddb3 100644
--- a/src/include/particle_simple.hxx
+++ b/src/include/particle_simple.hxx
@@ -50,13 +50,13 @@ public:
   template <typename FUNC>
   static void run(FUNC&& func)
   {
-    func("x", [](Particle& prt) { return &prt.x[0]; });
-    func("y", [](Particle& prt) { return &prt.x[1]; });
-    func("z", [](Particle& prt) { return &prt.x[2]; });
-    func("ux", [](Particle& prt) { return &prt.u[0]; });
-    func("uy", [](Particle& prt) { return &prt.u[1]; });
-    func("uz", [](Particle& prt) { return &prt.u[2]; });
-    func("kind", [](Particle& prt) { return &prt.kind; });
-    func("qni_wni", [](Particle& prt) { return &prt.qni_wni; });
+    func("x", [](auto& prt) { return &prt.x[0]; });
+    func("y", [](auto& prt) { return &prt.x[1]; });
+    func("z", [](auto& prt) { return &prt.x[2]; });
+    func("ux", [](auto& prt) { return &prt.u[0]; });
+    func("uy", [](auto& prt) { return &prt.u[1]; });
+    func("uz", [](auto& prt) { return &prt.u[2]; });
+    func("kind", [](auto& prt) { return &prt.kind; });
+    func("qni_wni", [](auto& prt) { return &prt.qni_wni; });
   }
 };
diff --git a/src/include/particles_simple.hxx b/src/include/particles_simple.hxx
index 5a68e6789e..0fbd773270 100644
--- a/src/include/particles_simple.hxx
+++ b/src/include/particles_simple.hxx
@@ -43,6 +43,7 @@ struct MparticlesStorage
   using PatchBuffer = std::vector<Particle>;
   using Buffers = std::vector<PatchBuffer>;
   using Range = Span<Particle>;
+  using ConstRange = Span<const Particle>;
   using iterator = typename Range::iterator;
   using const_iterator = typename Range::const_iterator;
 
@@ -91,10 +92,22 @@ struct MparticlesStorage
   }
 
   Range operator[](int p) { return {bufs_[p].data(), bufs_[p].size()}; }
+
+  ConstRange operator[](int p) const
+  {
+    return {bufs_[p].data(), bufs_[p].size()};
+  }
+
   Particle& at(int p, int n)
   {
     return bufs_[p][n];
   } // FIXME, ugly and not great for effciency
+
+  const Particle& at(int p, int n) const
+  {
+    return bufs_[p][n];
+  } // FIXME, ugly and not great for effciency
+
   void push_back(int p, const Particle& prt) { bufs_[p].push_back(prt); }
 
   Buffers& bndBuffers() { return bufs_; }
@@ -119,65 +132,7 @@ struct MparticlesSimple : MparticlesBase
   using BndBuffer = typename Storage::PatchBuffer;
   using BndBuffers = typename Storage::Buffers;
 
-  struct Patch
-  {
-    using iterator = typename Storage::iterator;
-    using const_iterator = typename Storage::const_iterator;
-
-    Patch(MparticlesSimple& mprts, int p) : mprts_(mprts), p_(p) {}
-
-    Patch(const Patch&) = delete;
-    Patch(Patch&&) = default;
-
-    Particle& operator[](int n) { return mprts_.storage_.at(p_, n); }
-    const Particle& operator[](int n) const
-    {
-      return mprts_.storage_.at(p_, n);
-    }
-
-    iterator begin() { return mprts_.storage_[p_].begin(); }
-    iterator end() { return mprts_.storage_[p_].end(); }
-    unsigned int size() const { return mprts_.storage_[p_].size(); }
-
-    void push_back(const Particle& new_prt)
-    {
-      // need to copy because we modify it
-      auto prt = new_prt;
-      checkInPatchMod(prt);
-      validCellIndex(prt);
-      mprts_.storage_.push_back(p_, prt);
-    }
-
-    void check() const
-    {
-      for (auto& prt : mprts_.storage_[p_]) {
-        mprts_.pi_.validCellIndex(prt.x());
-      }
-    }
-
-    // ParticleIndexer functionality
-    int cellPosition(real_t xi, int d) const
-    {
-      return mprts_.pi_.cellPosition(xi, d);
-    }
-    int validCellIndex(const Particle& prt) const
-    {
-      return mprts_.pi_.validCellIndex(prt.x);
-    }
-
-    void checkInPatchMod(Particle& prt) const
-    {
-      return mprts_.pi_.checkInPatchMod(prt.x);
-    }
-
-    const Grid_t& grid() const { return mprts_.grid(); }
-    const MparticlesSimple& mprts() const { return mprts_; }
-    int p() const { return p_; }
-
-  private:
-    MparticlesSimple& mprts_;
-    int p_;
-  };
+  using iterator = typename Storage::iterator;
 
   explicit MparticlesSimple(const Grid_t& grid)
     : MparticlesBase(grid),
@@ -197,11 +152,6 @@ struct MparticlesSimple : MparticlesBase
     storage_.reset(grid);
   }
 
-  Patch operator[](int p) const
-  {
-    return {const_cast<MparticlesSimple&>(*this), p};
-  } // FIXME, isn't actually const
-
   void reserve_all(const std::vector<uint>& n_prts_by_patch)
   {
     storage_.reserve_all(n_prts_by_patch);
@@ -219,6 +169,22 @@ struct MparticlesSimple : MparticlesBase
 
   const ParticleIndexer<real_t>& particleIndexer() const { return pi_; }
 
+  int cellPosition(int p, int n, int d) const
+  {
+    return pi_.cellPosition(this->at(p, n).x[d], d);
+  }
+
+  int validCellIndex(const Real3& x) const { return pi_.validCellIndex(x); }
+
+  void push_back(int p, const Particle& new_prt)
+  {
+    // need to copy because we modify it
+    auto prt = new_prt;
+    pi_.checkInPatchMod(prt.x);
+    validCellIndex(prt.x);
+    storage_.push_back(p, prt);
+  }
+
   InjectorSimple<MparticlesSimple> injector() { return {*this}; }
   ConstAccessor accessor() const
   {
@@ -228,13 +194,6 @@ struct MparticlesSimple : MparticlesBase
 
   BndBuffers& bndBuffers() { return storage_.bndBuffers(); }
 
-  void check() const
-  {
-    for (int p = 0; p < n_patches(); p++) {
-      (*this)[p].check();
-    }
-  }
-
   void dump(const std::string& filename)
   {
     FILE* file = fopen(filename.c_str(), "w");
@@ -254,6 +213,16 @@ struct MparticlesSimple : MparticlesBase
     fclose(file);
   }
 
+  Particle& at(int p, int n) { return storage_.at(p, n); }
+
+  const Particle& at(int p, int n) const { return storage_.at(p, n); }
+
+  iterator begin(int p) { return storage_[p].begin(); }
+
+  iterator end(int p) { return storage_[p].end(); }
+
+  unsigned int size(int p) const { return storage_[p].size(); }
+
   real_t prt_q(const Particle& prt) const { return grid().kinds[prt.kind].q; }
 
   real_t prt_m(const Particle& prt) const { return grid().kinds[prt.kind].m; }
diff --git a/src/include/particles_simple.inl b/src/include/particles_simple.inl
index eacd18066d..094d9c5556 100644
--- a/src/include/particles_simple.inl
+++ b/src/include/particles_simple.inl
@@ -53,13 +53,13 @@ public:
   template <typename FUNC>
   void operator()(const std::string& name, FUNC&& func)
   {
-    using Ret = typename std::remove_pointer<decltype(func(mprts_[0][0]))>::type;
+    using Ret = typename std::remove_const_t<
+      typename std::remove_pointer_t<decltype(func(mprts_.at(0, 0)))>>;
     std::vector<Ret> vec(mprts_.size());
     auto it = vec.begin();
     for (int p = 0; p < mprts_.n_patches(); p++) {
-      auto prts = mprts_[p];
-      for (int n = 0; n < prts.size(); n++) {
-        *it++ = *func(prts[n]);
+      for (int n = 0; n < mprts_.size(p); n++) {
+        *it++ = *func(mprts_.at(p, n));
       }
     }
 
@@ -83,15 +83,15 @@ public:
   template <typename FUNC>
   void operator()(const std::string& name, FUNC&& func)
   {
-    using Ret = typename std::remove_pointer<decltype(func(mprts_[0][0]))>::type;
+    using Ret =
+      typename std::remove_pointer<decltype(func(mprts_.at(0, 0)))>::type;
     std::vector<Ret> vec(mprts_.size());
     reader_.get<VariableByParticle>(name, vec, mprts_.grid(),
                                     kg::io::Mode::Blocking);
     auto it = vec.begin();
     for (int p = 0; p < mprts_.n_patches(); p++) {
-      auto prts = mprts_[p];
-      for (int n = 0; n < prts.size(); n++) {
-        *func(prts[n]) = *it++;
+      for (int n = 0; n < mprts_.size(p); n++) {
+        *func(mprts_.at(p, n)) = *it++;
       }
     }
   }
diff --git a/src/libpsc/mparticles.cxx b/src/libpsc/mparticles.cxx
index d0386be8c1..8076521e20 100644
--- a/src/libpsc/mparticles.cxx
+++ b/src/libpsc/mparticles.cxx
@@ -49,7 +49,6 @@ void psc_mparticles_check(MparticlesBase& mprts_base)
 
   for (int p = 0; p < grid.n_patches(); p++) {
     auto& patch = grid.patches[p];
-    auto&& prts = mprts[p];
 
     double xb[3], xe[3];
 
@@ -61,12 +60,13 @@ void psc_mparticles_check(MparticlesBase& mprts_base)
       xe[d] = patch.xb[d] + grid.ldims[d] * grid.domain.dx[d];
     }
 
-    for (auto prt : prts) {
-      if (prt.x[0] < 0.f || prt.x[0] >= xe[0] - xb[0] || // FIXME xz only!
-          prt.x[2] < 0.f || prt.x[2] >= xe[2] - xb[2]) {
+    for (auto prt_iter = mprts.begin(p); prt_iter != mprts.end(p); prt_iter++) {
+      if (prt_iter->x[0] < 0.f ||
+          prt_iter->x[0] >= xe[0] - xb[0] || // FIXME xz only!
+          prt_iter->x[2] < 0.f || prt_iter->x[2] >= xe[2] - xb[2]) {
         if (fail_cnt++ < 10) {
-          mprintf("FAIL: xi %g [%g:%g]\n", prt.x[0], 0., xe[0] - xb[0]);
-          mprintf("      zi %g [%g:%g]\n", prt.x[2], 0., xe[2] - xb[2]);
+          mprintf("FAIL: xi %g [%g:%g]\n", prt_iter->x[0], 0., xe[0] - xb[0]);
+          mprintf("      zi %g [%g:%g]\n", prt_iter->x[2], 0., xe[2] - xb[2]);
         }
       }
     }
diff --git a/src/libpsc/psc_balance/psc_balance_impl.hxx b/src/libpsc/psc_balance/psc_balance_impl.hxx
index 97084c4643..78164a4457 100644
--- a/src/libpsc/psc_balance/psc_balance_impl.hxx
+++ b/src/libpsc/psc_balance/psc_balance_impl.hxx
@@ -608,9 +608,8 @@ public:
 
       for (int pi = 0; pi < recv->nr_patches; pi++) {
         int p = recv->pi_to_patch[pi];
-        auto&& prts_new = mp_new[p];
-        int nn = prts_new.size() * (sizeof(Particle) / sizeof(real_t));
-        MPI_Irecv(&*prts_new.begin(), nn, mpi_dtype, recv->rank, pi, comm_,
+        int nn = mp_new.size(p) * (sizeof(Particle) / sizeof(real_t));
+        MPI_Irecv(&*mp_new.begin(p), nn, mpi_dtype, recv->rank, pi, comm_,
                   &recv_reqs[nr_recv_reqs++]);
       }
     }
@@ -629,10 +628,9 @@ public:
 
       for (int pi = 0; pi < send->nr_patches; pi++) {
         int p = send->pi_to_patch[pi];
-        auto&& prts_old = mp_old[p];
-        int nn = prts_old.size() * (sizeof(Particle) / sizeof(real_t));
+        int nn = mp_old.size(p) * (sizeof(Particle) / sizeof(real_t));
         // mprintf("A send -> %d tag %d (patch %d)\n", send->rank, pi, p);
-        MPI_Isend(&*prts_old.begin(), nn, mpi_dtype, send->rank, pi, comm_,
+        MPI_Isend(&*mp_old.begin(p), nn, mpi_dtype, send->rank, pi, comm_,
                   &send_reqs[nr_send_reqs++]);
       }
     }
@@ -646,12 +644,10 @@ public:
         continue;
       }
 
-      auto&& prts_old = mp_old[recv_info_[p].patch];
-      auto&& prts_new = mp_new[p];
-      assert(prts_old.size() == prts_new.size());
+      assert(mp_old.size(recv_info_[p].patch) == mp_new.size(p));
 #if 1
-      for (int n = 0; n < prts_new.size(); n++) {
-        prts_new[n] = prts_old[n];
+      for (int n = 0; n < mp_new.size(p); n++) {
+        mp_new.at(p, n) = mp_old.at(recv_info_[p].patch, n);
       }
 #else
       // FIXME, this needs at least a proper interface -- if not separately
diff --git a/src/libpsc/psc_collision/psc_collision_impl.hxx b/src/libpsc/psc_collision/psc_collision_impl.hxx
index eadf79688e..89c9887f66 100644
--- a/src/libpsc/psc_collision/psc_collision_impl.hxx
+++ b/src/libpsc/psc_collision/psc_collision_impl.hxx
@@ -23,7 +23,6 @@ struct CollisionHost
   using real_t = typename Mparticles::real_t;
   using MfieldsState = _MfieldsState;
   using Mfields = _Mfields;
-  using Particles = typename Mparticles::Patch;
   using Particle = typename Mparticles::Particle;
 
   enum
@@ -59,28 +58,26 @@ struct CollisionHost
     auto& grid = mprts.grid();
 
     for (int p = 0; p < mprts.n_patches(); p++) {
-      auto prts = mprts[p];
-
       const int* ldims = grid.ldims;
       int nr_cells = ldims[0] * ldims[1] * ldims[2];
       int* offsets = (int*)calloc(nr_cells + 1, sizeof(*offsets));
       struct psc_collision_stats stats_total = {};
 
-      find_cell_offsets(prts, offsets);
+      find_cell_offsets(mprts, p, offsets);
 
       auto F = make_Fields3d<dim_xyz>(mflds_stats_[p]);
       grid.Foreach_3d(0, 0, [&](int ix, int iy, int iz) {
         int c = (iz * ldims[1] + iy) * ldims[0] + ix;
 
-        update_rei_before(prts, offsets[c], offsets[c + 1], p, ix, iy, iz);
+        update_rei_before(mprts, p, offsets[c], offsets[c + 1], ix, iy, iz);
 
         struct psc_collision_stats stats = {};
         // mprintf("p %d ijk %d:%d:%d # %d\n", p, ix, iy, iz, offsets[c+1] -
         // offsets[c]);
         auto permute = randomize_in_cell(offsets[c], offsets[c + 1]);
-        collide_in_cell(prts, permute, &stats);
+        collide_in_cell(mprts, p, permute, &stats);
 
-        update_rei_after(prts, offsets[c], offsets[c + 1], p, ix, iy, iz);
+        update_rei_after(mprts, p, offsets[c], offsets[c + 1], ix, iy, iz);
 
         for (int s = 0; s < NR_STATS; s++) {
           F(s, ix, iy, iz) = stats.s[s];
@@ -144,14 +141,14 @@ struct CollisionHost
   // ----------------------------------------------------------------------
   // find_cell_offsets
 
-  static void find_cell_offsets(const Particles& prts, int offsets[])
+  static void find_cell_offsets(const Mparticles& mprts, int p, int offsets[])
   {
-    const int* ldims = prts.grid().ldims;
+    const int* ldims = mprts.grid().ldims;
     int last = 0;
     offsets[last] = 0;
-    int n_prts = prts.size();
+    int n_prts = mprts.size(p);
     for (int n = 0; n < n_prts; n++) {
-      int cell_index = prts.validCellIndex(prts[n]);
+      int cell_index = mprts.validCellIndex(mprts.at(p, n).x);
       assert(cell_index >= last);
       while (last < cell_index) {
         offsets[++last] = n;
@@ -176,17 +173,16 @@ struct CollisionHost
   // ----------------------------------------------------------------------
   // update_rei_before
 
-  void update_rei_before(const Particles& prts, int n_start, int n_end, int p,
+  void update_rei_before(const Mparticles& mprts, int p, int n_start, int n_end,
                          int i, int j, int k)
   {
-    real_t fnqs = prts.grid().norm.fnqs;
+    real_t fnqs = mprts.grid().norm.fnqs;
     auto F = make_Fields3d<dim_xyz>(mflds_rei_[p]);
     F(0, i, j, k) = 0.;
     F(1, i, j, k) = 0.;
     F(2, i, j, k) = 0.;
-    auto& mprts = prts.mprts();
     for (int n = n_start; n < n_end; n++) {
-      const auto& prt = prts[n];
+      const auto& prt = mprts.at(p, n);
       auto fac = mprts.prt_m(prt) * mprts.prt_w(prt) * fnqs;
       F(0, i, j, k) -= prt.u[0] * fac;
       F(1, i, j, k) -= prt.u[1] * fac;
@@ -197,14 +193,14 @@ struct CollisionHost
   // ----------------------------------------------------------------------
   // update_rei_after
 
-  void update_rei_after(const Particles& prts, int n_start, int n_end, int p,
+  void update_rei_after(const Mparticles& mprts, int p, int n_start, int n_end,
                         int i, int j, int k)
   {
-    real_t fnqs = prts.grid().norm.fnqs, dt = prts.grid().dt;
+    real_t fnqs = mprts.grid().norm.fnqs;
+    real_t dt = mprts.grid().dt;
     auto F = make_Fields3d<dim_xyz>(mflds_rei_[p]);
-    auto& mprts = prts.mprts();
     for (int n = n_start; n < n_end; n++) {
-      const auto& prt = prts[n];
+      const auto& prt = mprts.at(p, n);
       auto fac = mprts.prt_m(prt) * mprts.prt_w(prt) * fnqs;
       F(0, i, j, k) += prt.u[0] * fac;
       F(1, i, j, k) += prt.u[1] * fac;
@@ -218,19 +214,19 @@ struct CollisionHost
   // ----------------------------------------------------------------------
   // collide_in_cell
 
-  void collide_in_cell(Particles& prts, const std::vector<int>& permute,
+  void collide_in_cell(Mparticles& mprts, int p,
+                       const std::vector<int>& permute,
                        struct psc_collision_stats* stats)
   {
-    const auto& grid = prts.grid();
+    const auto& grid = mprts.grid();
     int nn = permute.size();
 
     if (nn < 2) { // can't collide only one (or zero) particles
       return;
     }
 
-    const auto& mprts = prts.mprts();
     // all particles need to have same weight!
-    real_t wni = mprts.prt_w(prts[permute[0]]);
+    real_t wni = mprts.prt_w(mprts.at(p, permute[0]));
     real_t nudt1 = wni * grid.norm.cori * nn * this->interval_ * grid.dt * nu_;
 
     real_t* nudts = (real_t*)malloc((nn / 2 + 2) * sizeof(*nudts));
@@ -238,25 +234,24 @@ struct CollisionHost
 
     int n = 0;
     if (nn % 2 == 1) { // odd # of particles: do 3-collision
-      nudts[cnt++] = do_bc(prts, permute[0], permute[1], .5 * nudt1);
-      nudts[cnt++] = do_bc(prts, permute[0], permute[2], .5 * nudt1);
-      nudts[cnt++] = do_bc(prts, permute[1], permute[2], .5 * nudt1);
+      nudts[cnt++] = do_bc(mprts, p, permute[0], permute[1], .5 * nudt1);
+      nudts[cnt++] = do_bc(mprts, p, permute[0], permute[2], .5 * nudt1);
+      nudts[cnt++] = do_bc(mprts, p, permute[1], permute[2], .5 * nudt1);
       n = 3;
     }
     for (; n < nn; n += 2) { // do remaining particles as pair
-      nudts[cnt++] = do_bc(prts, permute[n], permute[n + 1], nudt1);
+      nudts[cnt++] = do_bc(mprts, p, permute[n], permute[n + 1], nudt1);
     }
 
     calc_stats(stats, nudts, cnt);
     free(nudts);
   }
 
-  real_t do_bc(Particles& prts, int n1, int n2, real_t nudt1)
+  real_t do_bc(Mparticles& mprts, int p, int n1, int n2, real_t nudt1)
   {
     Rng rng;
-    const auto& mprts = prts.mprts();
     BinaryCollision<Mparticles, Particle> bc(mprts);
-    return bc(prts[n1], prts[n2], nudt1, rng);
+    return bc(mprts.at(p, n1), mprts.at(p, n2), nudt1, rng);
   }
 
   int interval() const { return interval_; }
diff --git a/src/libpsc/psc_heating/psc_heating_impl.hxx b/src/libpsc/psc_heating/psc_heating_impl.hxx
index f82c1a1142..07c794fe81 100644
--- a/src/libpsc/psc_heating/psc_heating_impl.hxx
+++ b/src/libpsc/psc_heating/psc_heating_impl.hxx
@@ -55,18 +55,17 @@ struct Heating__ : HeatingBase
   void operator()(Mparticles& mprts)
   {
     for (int p = 0; p < mprts.n_patches(); p++) {
-      auto&& prts = mprts[p];
       auto& patch = mprts.grid().patches[p];
-      for (auto& prt : prts) {
-
+      for (auto prt_iter = mprts.begin(p); prt_iter != mprts.end(p);
+           prt_iter++) {
         double xx[3] = {
-          prt.x[0] + patch.xb[0],
-          prt.x[1] + patch.xb[1],
-          prt.x[2] + patch.xb[2],
+          prt_iter->x[0] + patch.xb[0],
+          prt_iter->x[1] + patch.xb[1],
+          prt_iter->x[2] + patch.xb[2],
         };
-        double H = get_H_(xx, prt.kind);
+        double H = get_H_(xx, prt_iter->kind);
         if (H > 0.f) {
-          kick_particle(prt, H);
+          kick_particle(*prt_iter, H);
         }
       }
     }
diff --git a/src/libpsc/psc_output_particles/output_particles_hdf5_impl.hxx b/src/libpsc/psc_output_particles/output_particles_hdf5_impl.hxx
index af1798b7af..1417bfd676 100644
--- a/src/libpsc/psc_output_particles/output_particles_hdf5_impl.hxx
+++ b/src/libpsc/psc_output_particles/output_particles_hdf5_impl.hxx
@@ -433,7 +433,6 @@ struct OutputParticlesHdf5
 {
   using writer_type = OutputParticlesWriterHDF5;
   using writer_particles_type = writer_type::particles_type;
-  using Particles = typename Mparticles::Patch;
 
   OutputParticlesHdf5(const Grid_t& grid, const OutputParticlesParams& params)
     : lo_{params.lo},
@@ -454,15 +453,15 @@ struct OutputParticlesHdf5
            kind;
   };
 
-  static inline int get_sort_index(Particles& prts, int n)
+  static inline int get_sort_index(Mparticles& mprts, int p, int n)
   {
-    const Grid_t& grid = prts.grid();
+    const Grid_t& grid = mprts.grid();
     const int* ldims = grid.ldims;
-    const auto& prt = prts[n];
+    const auto& prt = mprts.at(p, n);
 
     Int3 pos;
     for (int d = 0; d < 3; d++) {
-      pos[d] = prts.cellPosition(prt.x[d], d);
+      pos[d] = mprts.cellPosition(p, n, d);
       // FIXME, this is hoping that reason is that we were just on the right
       // bnd...
       if (pos[d] == ldims[d])
@@ -488,16 +487,15 @@ struct OutputParticlesHdf5
       const int* ldims = mprts.grid().ldims;
       int nr_indices = ldims[0] * ldims[1] * ldims[2] * nr_kinds;
       off[p].resize(nr_indices + 1);
-      auto&& prts = mprts[p];
-      unsigned int n_prts = prts.size();
+      unsigned int n_prts = mprts.size(p);
       std::vector<int> particle_indices;
       particle_indices.reserve(n_prts);
 
       // counting sort to get map
       for (int n = 0; n < n_prts; n++) {
-        if (selector(prts[n])) {
+        if (selector(mprts.at(p, n))) {
           particle_indices.push_back(n);
-          int si = get_sort_index(prts, n);
+          int si = get_sort_index(mprts, p, n);
           off[p][si]++;
         }
       }
@@ -514,7 +512,7 @@ struct OutputParticlesHdf5
       // sort a map only, not the actual particles
       map[p].resize(n_prts);
       for (auto n : particle_indices) {
-        int si = get_sort_index(prts, n);
+        int si = get_sort_index(mprts, p, n);
         map[p][off2[si]++] = n;
       }
     }
diff --git a/src/libpsc/psc_particles_impl.cxx b/src/libpsc/psc_particles_impl.cxx
index 800a7d51bf..74389d4e99 100644
--- a/src/libpsc/psc_particles_impl.cxx
+++ b/src/libpsc/psc_particles_impl.cxx
@@ -30,11 +30,9 @@ void psc_mparticles_copy(MP_FROM& mp_from, MP_TO& mp_to)
   mp_to.clear();
 
   for (int p = 0; p < mp_to.n_patches(); p++) {
-    auto&& prts_from = mp_from[p];
-    auto&& prts_to = mp_to[p];
-    int n_prts = prts_from.size();
+    int n_prts = mp_from.size(p);
     for (int n = 0; n < n_prts; n++) {
-      prts_to.push_back(convert(prts_from[n], mp_from.grid()));
+      mp_to.push_back(p, convert(mp_from.at(p, n), mp_from.grid()));
     }
   }
 }
diff --git a/src/libpsc/psc_sort/psc_sort_impl.hxx b/src/libpsc/psc_sort/psc_sort_impl.hxx
index 9c64412dee..ee985d746f 100644
--- a/src/libpsc/psc_sort/psc_sort_impl.hxx
+++ b/src/libpsc/psc_sort/psc_sort_impl.hxx
@@ -71,16 +71,15 @@ struct SortCountsort2
   void operator()(Mparticles& mprts)
   {
     for (int p = 0; p < mprts.n_patches(); p++) {
-      auto&& prts = mprts[p];
-      unsigned int n_prts = prts.size();
+      unsigned int n_prts = mprts.size(p);
 
       unsigned int n_cells = mprts.pi_.n_cells_;
       unsigned int* cnis = new unsigned int[n_prts];
       // FIXME, might as well merge counting here, too
       int i = 0;
-      for (auto prt_iter = prts.begin(); prt_iter != prts.end();
+      for (auto prt_iter = mprts.begin(p); prt_iter != mprts.end(p);
            ++prt_iter, ++i) {
-        cnis[i] = prts.validCellIndex(*prt_iter);
+        cnis[i] = mprts.validCellIndex(prt_iter->x);
       }
 
       unsigned int* cnts = new unsigned int[n_cells]{};
@@ -108,13 +107,14 @@ struct SortCountsort2
         while (i + n < n_prts && cnis[i + n] == cni) {
           n++;
         }
-        memcpy(&particles2[cnts[cni]], &prts[i], n * sizeof(*particles2));
+        memcpy(&particles2[cnts[cni]], &mprts.at(p, i),
+               n * sizeof(*particles2));
         cnts[cni] += n;
         i += n - 1;
       }
 
       // back to in-place
-      memcpy(&*prts.begin(), particles2, n_prts * sizeof(*particles2));
+      memcpy(&*mprts.begin(p), particles2, n_prts * sizeof(*particles2));
 
       delete[] particles2;
       delete[] cnis;
diff --git a/src/libpsc/tests/test_mparticles.cxx b/src/libpsc/tests/test_mparticles.cxx
index cc8e362e22..770418669f 100644
--- a/src/libpsc/tests/test_mparticles.cxx
+++ b/src/libpsc/tests/test_mparticles.cxx
@@ -503,9 +503,10 @@ TEST(TestSetupParticles, Id)
   EXPECT_EQ(mprts.size(), n_cells * kinds.size() * prm.nicell);
 
   psc::particle::Id cnt = 0;
+
   for (int p = 0; p < mprts.n_patches(); p++) {
-    for (auto& prt : mprts[p]) {
-      EXPECT_EQ(prt.id(), cnt++);
+    for (auto prt_iter = mprts.begin(p); prt_iter != mprts.end(p); prt_iter++) {
+      EXPECT_EQ(prt_iter->id(), cnt++);
     }
   }
 }
@@ -532,8 +533,8 @@ TEST(TestSetupParticles, Tag)
   EXPECT_EQ(mprts.size(), n_cells * kinds.size() * prm.nicell);
 
   for (int p = 0; p < mprts.n_patches(); p++) {
-    for (auto& prt : mprts[p]) {
-      EXPECT_EQ(prt.tag(), prt.kind * 10);
+    for (auto prt_iter = mprts.begin(p); prt_iter != mprts.end(p); prt_iter++) {
+      EXPECT_EQ(prt_iter->tag(), prt_iter->kind * 10);
     }
   }
 }
diff --git a/src/libpsc/vpic/psc_particles_vpic.cxx b/src/libpsc/vpic/psc_particles_vpic.cxx
index 4022751173..ddca96e932 100644
--- a/src/libpsc/vpic/psc_particles_vpic.cxx
+++ b/src/libpsc/vpic/psc_particles_vpic.cxx
@@ -52,8 +52,8 @@ void copy_to(MparticlesBase& mprts_from_base, MparticlesBase& mprts_to_base)
         auto kind = sp->id;
         auto qni_wni =
           float(vprt.w * dVi) * float(mprts_to.grid().kinds[kind].q);
-        mprts_to[p].push_back(
-          {x, u, qni_wni, kind, psc::particle::Id{}, psc::particle::Tag{}});
+        mprts_to.push_back(
+          p, {x, u, qni_wni, kind, psc::particle::Id{}, psc::particle::Tag{}});
       }
     }
   }
diff --git a/src/psc_flatfoil_yz.cxx b/src/psc_flatfoil_yz.cxx
index ff2910e74f..f2de74e21a 100644
--- a/src/psc_flatfoil_yz.cxx
+++ b/src/psc_flatfoil_yz.cxx
@@ -436,12 +436,11 @@ void run()
       // fix up renumbered particle kind
       auto& mp = mprts.template get_as<MparticlesSingle>();
       for (int p = 0; p < mp.n_patches(); p++) {
-        auto prts = mp[p];
-        for (int n = 0; n < prts.size(); n++) {
-          if (prts[n].kind == 0) {
-            prts[n].kind = MY_ION;
-          } else if (prts[n].kind == 1) {
-            prts[n].kind = MY_ELECTRON;
+        for (int n = 0; n < mprts.size(p); n++) {
+          if (mprts.at(p, n).kind == 0) {
+            mprts.at(p, n).kind = MY_ION;
+          } else if (mprts.at(p, n).kind == 1) {
+            mprts.at(p, n).kind = MY_ELECTRON;
           } else {
             assert(0);
           }