From 1227c91cfc12c3a7298966cd228429d1c1bd770e Mon Sep 17 00:00:00 2001
From: Axel Huebl <axel.huebl@plasma.ninja>
Date: Tue, 30 Jan 2024 14:01:04 -0800
Subject: [PATCH] Particle Container to Pure SoA Again

Transition to new, purely SoA particle containers.

This was originally merged in #3850 and reverted in #4652, since
we discovered issues loosing particles & laser particles on GPU.
---
 Docs/source/developers/amrex_basics.rst       |   2 +-
 Docs/source/developers/dimensionality.rst     |   4 +-
 Docs/source/developers/particles.rst          |  18 +-
 Docs/source/usage/workflows/python_extend.rst |  28 +-
 .../particle_data_python/PICMI_inputs_2d.py   |   6 +-
 .../PICMI_inputs_prev_pos_2d.py               |   6 +-
 .../PICMI_inputs_runtime_component_analyze.py |   6 +-
 Python/pywarpx/_libwarpx.py                   |   6 +-
 Python/pywarpx/particle_containers.py         | 255 ++++++++----------
 .../BackTransformParticleFunctor.H            |  16 +-
 .../FlushFormats/FlushFormatAscent.cpp        |   3 -
 .../FlushFormats/FlushFormatCheckpoint.cpp    |   9 +-
 .../FlushFormats/FlushFormatPlotfile.cpp      |  13 +-
 Source/Diagnostics/ParticleIO.cpp             |   2 +-
 .../Diagnostics/ReducedDiags/FieldProbe.cpp   |   6 +-
 .../FieldProbeParticleContainer.H             |  20 +-
 .../FieldProbeParticleContainer.cpp           |  36 +--
 .../ReducedDiags/LoadBalanceCosts.cpp         |   3 +-
 Source/Diagnostics/WarpXOpenPMD.H             |   4 +-
 Source/Diagnostics/WarpXOpenPMD.cpp           | 139 ++--------
 .../ParticleBoundaryProcess.H                 |   5 +-
 Source/EmbeddedBoundary/ParticleScraper.H     |   3 +-
 .../BinaryCollision/BinaryCollision.H         |  16 +-
 .../Coulomb/PairWiseCoulombCollisionFunc.H    |   7 +-
 .../Collision/BinaryCollision/DSMC/DSMC.H     |   3 +-
 .../DSMC/SplitAndScatterFunc.H                |  30 ++-
 .../NuclearFusion/NuclearFusionFunc.H         |  14 +-
 .../ProtonBoronFusionInitializeMomentum.H     |  10 +-
 .../TwoProductFusionInitializeMomentum.H      |   4 +-
 .../BinaryCollision/ParticleCreationFunc.H    |  60 +++--
 .../BinaryCollision/ShuffleFisherYates.H      |   2 +-
 .../Particles/Deposition/ChargeDeposition.H   |   2 +-
 .../Particles/Deposition/CurrentDeposition.H  |   2 +-
 .../ElementaryProcess/QEDPairGeneration.H     |   2 +-
 .../ElementaryProcess/QEDPhotonEmission.H     |  16 +-
 Source/Particles/LaserParticleContainer.H     |   7 +-
 .../NamedComponentParticleContainer.H         |  52 ++--
 Source/Particles/ParticleBoundaryBuffer.cpp   |   4 +-
 Source/Particles/ParticleCreation/SmartCopy.H |   5 +-
 .../Particles/ParticleCreation/SmartCreate.H  |  22 +-
 .../Particles/ParticleCreation/SmartUtils.H   |   9 +-
 Source/Particles/PhysicalParticleContainer.H  |   8 +-
 .../Particles/PhysicalParticleContainer.cpp   | 129 ++++-----
 Source/Particles/Pusher/GetAndSetPosition.H   | 152 +++++++----
 .../Particles/Resampling/LevelingThinning.cpp |   5 +-
 Source/Particles/Sorting/Partition.cpp        |   4 +-
 Source/Particles/Sorting/SortingUtils.H       |  41 ++-
 Source/Particles/Sorting/SortingUtils.cpp     |   2 +-
 Source/Particles/WarpXParticleContainer.H     |  27 +-
 Source/Particles/WarpXParticleContainer.cpp   |  89 +++---
 .../Particles/ParticleBoundaryBuffer.cpp      |  10 +-
 .../PinnedMemoryParticleContainer.cpp         |   2 +-
 .../Particles/WarpXParticleContainer.cpp      |  18 +-
 Source/Utils/ParticleUtils.H                  |   7 +-
 Source/Utils/ParticleUtils.cpp                |  26 +-
 Source/ablastr/particles/IndexHandling.H      |  41 ---
 Source/ablastr/particles/ParticleMoments.H    |  25 +-
 57 files changed, 700 insertions(+), 743 deletions(-)
 delete mode 100644 Source/ablastr/particles/IndexHandling.H

diff --git a/Docs/source/developers/amrex_basics.rst b/Docs/source/developers/amrex_basics.rst
index 577a6547bb5..64ad71af06c 100644
--- a/Docs/source/developers/amrex_basics.rst
+++ b/Docs/source/developers/amrex_basics.rst
@@ -13,7 +13,7 @@ WarpX is built on the Adaptive Mesh Refinement (AMR) library `AMReX <https://git
 
 * ``amrex::MultiFab``: Collection of `FAB` (= ``FArrayBox``) on a single AMR level, distributed over MPI ranks. The concept of `ghost cells` is defined at the ``MultiFab`` level.
 
-* ``amrex::ParticleContainer``: A collection of particles, typically for particles of a physical species. Particles in a ``ParticleContainer`` are organized per ``Box``. Particles in a ``Box`` are organized per tile (this feature is off when running on GPU). Particles within a tile are stored in several structures, each being contiguous in memory: (i) an Array-Of-Struct (AoS) (often called `data`, they are the 3D position, the particle ID and the index of the CPU owning the particle), where the Struct is an ``amrex::Particle`` and (ii) Struct-Of-Arrays (SoA) for extra variables (often called ``attribs``, in WarpX they are the momentum, field on particle etc.).
+* ``amrex::ParticleContainer``: A collection of particles, typically for particles of a physical species. Particles in a ``ParticleContainer`` are organized per ``Box``. Particles in a ``Box`` are organized per tile (this feature is off when running on GPU). Particles within a tile are stored in several structures, each being contiguous in memory: (i) a Struct-of-Array (SoA) for ``amrex::ParticleReal`` data such as positions, weight, momentum, etc., (ii) a Struct-of-Array (SoA) for ``int`` data, such as ionization levels, and (iii) a Struct-of-Array (SoA) for a ``uint64_t`` unique identifier index per particle (containing a 40bit id and 24bit cpu sub-identifier as assigned at particle creation time). This id is also used to check if a particle is active/valid or marked for removal.
 
 The simulation domain is decomposed in several ``Box``, and each MPI rank owns (and performs operations on) the fields and particles defined on a few of these ``Box``, but has the metadata of all of them. For convenience, AMReX provides iterators, to easily iterate over all ``FArrayBox`` (or even tile-by-tile, optionally) in a ``MultiFab`` own by the MPI rank (``MFIter``), or over all particles in a ``ParticleContainer`` on a per-box basis (``ParIter``, or its derived class ``WarpXParIter``). These are respectively done in loops like:
 
diff --git a/Docs/source/developers/dimensionality.rst b/Docs/source/developers/dimensionality.rst
index e775d31a6df..52493e91e38 100644
--- a/Docs/source/developers/dimensionality.rst
+++ b/Docs/source/developers/dimensionality.rst
@@ -49,12 +49,12 @@ WarpX axis labels     ``x, y, z``   ``x, z``    ``z``        ``x, z``
 --------------------  -----------  -----------  -----------  -----------
 *Particles*
 ------------------------------------------------------------------------
-AMReX AoS ``.pos()``  ``0, 1, 2``  ``0, 1``     ``0``        ``0, 1``
+AMReX ``.pos()``      ``0, 1, 2``  ``0, 1``     ``0``        ``0, 1``
 WarpX position names  ``x, y, z``  ``x, z``     ``z``        ``r, z``
 extra SoA attribute                                          ``theta``
 ====================  ===========  ===========  ===========  ===========
 
-Please see the following sections for particle AoS and SoA details.
+Please see the following sections for particle SoA details.
 
 Conventions
 -----------
diff --git a/Docs/source/developers/particles.rst b/Docs/source/developers/particles.rst
index c1a711e0a01..a0182f3845f 100644
--- a/Docs/source/developers/particles.rst
+++ b/Docs/source/developers/particles.rst
@@ -109,24 +109,26 @@ Particle attributes
 -------------------
 
 WarpX adds the following particle attributes by default to WarpX particles.
-These attributes are either stored in an Array-of-Struct (AoS) or Struct-of-Array (SoA) location of the AMReX particle containers.
+These attributes are stored in Struct-of-Array (SoA) locations of the AMReX particle containers: one SoA for ``amrex::ParticleReal`` attributes, one SoA for ``int`` attributes and one SoA for a ``uint64_t`` global particle index per particle.
 The data structures for those are either pre-described at compile-time (CT) or runtime (RT).
 
-====================  ================  ==================================  ===== ==== =====================
+====================  ================  ==================================  ===== ==== ======================
 Attribute name        ``int``/``real``  Description                         Where When Notes
-====================  ================  ==================================  ===== ==== =====================
-``position_x/y/z``    ``real``          Particle position.                  AoS   CT
-``cpu``               ``int``           CPU index where the particle        AoS   CT
+====================  ================  ==================================  ===== ==== ======================
+``position_x/y/z``    ``real``          Particle position.                  SoA   CT
+``weight``            ``real``          Particle position.                  SoA   CT
+``momentum_x/y/z``    ``real``          Particle position.                  SoA   CT
+``id``                ``amrex::Long``   CPU-local particle index            SoA   CT   First 40 bytes of
+                                        where the particle was created.                idcpu
+``cpu``               ``int``           CPU index where the particle        SoA   CT   Last 24 bytes of idcpu
                                         was created.
-``id``                ``int``           CPU-local particle index            AoS   CT
-                                        where the particle was created.
 ``ionizationLevel``   ``int``           Ion ionization level                SoA   RT   Added when ionization
                                                                                        physics is used.
 ``opticalDepthQSR``   ``real``          QED: optical depth of the Quantum-  SoA   RT   Added when PICSAR QED
                                         Synchrotron process                            physics is used.
 ``opticalDepthBW``    ``real``          QED: optical depth of the Breit-    SoA   RT   Added when PICSAR QED
                                         Wheeler process                                physics is used.
-====================  ================  ==================================  ===== ==== =====================
+====================  ================  ==================================  ===== ==== ======================
 
 WarpX allows extra runtime attributes to be added to particle containers (through ``AddRealComp("attrname")`` or ``AddIntComp("attrname")``).
 The attribute name can then be used to access the values of that attribute.
diff --git a/Docs/source/usage/workflows/python_extend.rst b/Docs/source/usage/workflows/python_extend.rst
index 6c9286c02ce..d3049508da6 100644
--- a/Docs/source/usage/workflows/python_extend.rst
+++ b/Docs/source/usage/workflows/python_extend.rst
@@ -190,6 +190,7 @@ Particles
    @callfromafterstep
    def my_after_step_callback():
        warpx = sim.extension.warpx
+       Config = sim.extension.Config
 
        # data access
        multi_pc = warpx.multi_particle_container()
@@ -200,27 +201,27 @@ Particles
        for lvl in range(pc.finest_level + 1):
            # get every local chunk of particles
            for pti in pc.iterator(pc, level=lvl):
-               # default layout: AoS with positions and cpuid
-               aos = pti.aos().to_numpy()
-
-               # additional compile-time and runtime attributes in SoA format
-               soa = pti.soa().to_numpy()
+               # compile-time and runtime attributes in SoA format
+               soa = pti.soa().to_cupy() if Config.have_gpu else \
+                     pti.soa().to_numpy()
 
                # notes:
                # Only the next lines are the "HOT LOOP" of the computation.
-               # For efficiency, use numpy array operation for speed on CPUs.
-               # For GPUs use .to_cupy() above and compute with cupy or numba.
+               # For speed, use array operation.
 
                # write to all particles in the chunk
-               # note: careful, if you change particle positions, you need to
+               # note: careful, if you change particle positions, you might need to
                #       redistribute particles before continuing the simulation step
-               # aos[()]["x"] = 0.30
-               # aos[()]["y"] = 0.35
-               # aos[()]["z"] = 0.40
+               soa.real[0][()] = 0.30  # x
+               soa.real[1][()] = 0.35  # y
+               soa.real[2][()] = 0.40  # z
 
-               for soa_real in soa.real:
+               # all other attributes: weight, momentum x, y, z, ...
+               for soa_real in soa.real[3:]:
                    soa_real[()] = 42.0
 
+               # by default empty unless ionization or QED physics is used
+               # or other runtime attributes were added manually
                for soa_int in soa.int:
                    soa_int[()] = 12
 
@@ -252,7 +253,8 @@ Particles can be added to the simulation at specific positions and with specific
 .. autoclass:: pywarpx.particle_containers.ParticleContainerWrapper
    :members:
 
-The ``get_particle_structs()`` and ``get_particle_arrays()`` functions are called
+The ``get_particle_real_arrays()``, ``get_particle_int_arrays()`` and
+``get_particle_idcpu_arrays()`` functions are called
 by several utility functions of the form ``get_particle_{comp_name}`` where
 ``comp_name`` is one of ``x``, ``y``, ``z``, ``r``, ``theta``, ``id``, ``cpu``,
 ``weight``, ``ux``, ``uy`` or ``uz``.
diff --git a/Examples/Tests/particle_data_python/PICMI_inputs_2d.py b/Examples/Tests/particle_data_python/PICMI_inputs_2d.py
index a4b7d9e134e..572871b8ed5 100755
--- a/Examples/Tests/particle_data_python/PICMI_inputs_2d.py
+++ b/Examples/Tests/particle_data_python/PICMI_inputs_2d.py
@@ -153,10 +153,10 @@ def add_particles():
 ##########################
 
 assert (elec_wrapper.nps == 270 / (2 - args.unique))
-assert (elec_wrapper.particle_container.get_comp_index('w') == 0)
-assert (elec_wrapper.particle_container.get_comp_index('newPid') == 4)
+assert (elec_wrapper.particle_container.get_comp_index('w') == 2)
+assert (elec_wrapper.particle_container.get_comp_index('newPid') == 6)
 
-new_pid_vals = elec_wrapper.get_particle_arrays('newPid', 0)
+new_pid_vals = elec_wrapper.get_particle_real_arrays('newPid', 0)
 for vals in new_pid_vals:
     assert np.allclose(vals, 5)
 
diff --git a/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py b/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py
index 1becd4464e7..5de9879f0f8 100755
--- a/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py
+++ b/Examples/Tests/particle_data_python/PICMI_inputs_prev_pos_2d.py
@@ -120,12 +120,12 @@
 elec_count = elec_wrapper.nps
 
 # check that the runtime attributes have the right indices
-assert (elec_wrapper.particle_container.get_comp_index('prev_x') == 4)
-assert (elec_wrapper.particle_container.get_comp_index('prev_z') == 5)
+assert (elec_wrapper.particle_container.get_comp_index('prev_x') == 6)
+assert (elec_wrapper.particle_container.get_comp_index('prev_z') == 7)
 
 # sanity check that the prev_z values are reasonable and
 # that the correct number of values are returned
-prev_z_vals = elec_wrapper.get_particle_arrays('prev_z', 0)
+prev_z_vals = elec_wrapper.get_particle_real_arrays('prev_z', 0)
 running_count = 0
 
 for z_vals in prev_z_vals:
diff --git a/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py b/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py
index 706dedb6959..32c9f4e5808 100755
--- a/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py
+++ b/Examples/Tests/restart/PICMI_inputs_runtime_component_analyze.py
@@ -158,10 +158,10 @@ def add_particles():
 ##########################
 
 assert electron_wrapper.nps == 90
-assert electron_wrapper.particle_container.get_comp_index("w") == 0
-assert electron_wrapper.particle_container.get_comp_index("newPid") == 4
+assert electron_wrapper.particle_container.get_comp_index("w") == 2
+assert electron_wrapper.particle_container.get_comp_index("newPid") == 6
 
-new_pid_vals = electron_wrapper.get_particle_arrays("newPid", 0)
+new_pid_vals = electron_wrapper.get_particle_real_arrays("newPid", 0)
 for vals in new_pid_vals:
     assert np.allclose(vals, 5)
 
diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py
index fa2044d4240..1eb410e65b8 100755
--- a/Python/pywarpx/_libwarpx.py
+++ b/Python/pywarpx/_libwarpx.py
@@ -5,8 +5,8 @@
 #
 # NOTE: We will reduce the libwarpx.py level of abstraction eventually!
 # Please add new functionality directly to pybind11-bound modules
-# and call them via sim.extension.libwarpx_so. ... and sim.extension.warpx.
-# ... from user code.
+# and call them via sim.extension.libwarpx_so. ... and sim.extension.Config and
+# sim.extension.warpx. ... from user code.
 #
 # Authors: Axel Huebl, Andrew Myers, David Grote, Remi Lehe, Weiqun Zhang
 #
@@ -108,6 +108,8 @@ def load_library(self):
                 from . import warpx_pybind_3d as cxx_3d
                 self.libwarpx_so = cxx_3d
                 self.dim = 3
+
+            self.Config = self.libwarpx_so.Config
         except ImportError:
             raise Exception(f"Dimensionality '{self.geometry_dim}' was not compiled in this Python install. Please recompile with -DWarpX_DIMS={_dims}")
 
diff --git a/Python/pywarpx/particle_containers.py b/Python/pywarpx/particle_containers.py
index 273a981f4bd..72a675ec43c 100644
--- a/Python/pywarpx/particle_containers.py
+++ b/Python/pywarpx/particle_containers.py
@@ -117,8 +117,10 @@ def add_particles(self, x=None, y=None, z=None, ux=None, uy=None,
                 kwargs[key] = np.full(maxlen, val)
 
         # --- The number of built in attributes
+        # --- The positions
+        built_in_attrs = libwarpx.dim
         # --- The three velocities
-        built_in_attrs = 3
+        built_in_attrs += 3
         if libwarpx.geometry_dim == 'rz':
             # --- With RZ, there is also theta
             built_in_attrs += 1
@@ -188,28 +190,25 @@ def add_real_comp(self, pid_name, comm=True):
         self.particle_container.add_real_comp(pid_name, comm)
 
 
-    def get_particle_structs(self, level, copy_to_host=False):
+    def get_particle_real_arrays(self, comp_name, level, copy_to_host=False):
         '''
-        This returns a list of numpy or cupy arrays containing the particle struct data
-        on each tile for this process. The particle data is represented as a structured
-        array and contains the particle 'x', 'y', 'z', and 'idcpu'.
+        This returns a list of numpy or cupy arrays containing the particle real array data
+        on each tile for this process.
 
-        Unless copy_to_host is specified, the data for the structs are
-        not copied, but share the underlying memory buffer with WarpX. The
+        Unless copy_to_host is specified, the data for the arrays are not
+        copied, but share the underlying memory buffer with WarpX. The
         arrays are fully writeable.
 
-        Note that cupy does not support structs:
-        https://github.com/cupy/cupy/issues/2031
-        and will return arrays of binary blobs for the AoS (DP: ``"|V24"``). If copied
-        to host via copy_to_host, we correct for the right numpy AoS type.
-
         Parameters
         ----------
 
-        level        : int
+        comp_name      : str
+            The component of the array data that will be returned
+
+        level          : int
             The refinement level to reference (default=0)
 
-        copy_to_host : bool
+        copy_to_host   : bool
             For GPU-enabled runs, one can either return the GPU
             arrays (the default) or force a device-to-host copy.
 
@@ -217,30 +216,29 @@ def get_particle_structs(self, level, copy_to_host=False):
         -------
 
         List of arrays
-            The requested particle struct data
+            The requested particle array data
         '''
-        particle_data = []
+        comp_idx = self.particle_container.get_comp_index(comp_name)
+
+        data_array = []
         for pti in libwarpx.libwarpx_so.WarpXParIter(self.particle_container, level):
+            soa = pti.soa()
+            idx = soa.GetRealData(comp_idx)
             if copy_to_host:
-                particle_data.append(pti.aos().to_numpy(copy=True))
+                data_array.append(idx.to_numpy(copy=True))
             else:
-                if libwarpx.amr.Config.have_gpu:
-                    libwarpx.amr.Print(
-                        "get_particle_structs: cupy does not yet support structs. "
-                        "https://github.com/cupy/cupy/issues/2031"
-                        "Did you mean copy_to_host=True?"
-                    )
                 xp, cupy_status = load_cupy()
                 if cupy_status is not None:
                     libwarpx.amr.Print(cupy_status)
-                aos_arr = xp.array(pti.aos(), copy=False)   # void blobs for cupy
-                particle_data.append(aos_arr)
-        return particle_data
+
+                data_array.append(xp.array(idx, copy=False))
+
+        return data_array
 
 
-    def get_particle_arrays(self, comp_name, level, copy_to_host=False):
+    def get_particle_int_arrays(self, comp_name, level, copy_to_host=False):
         '''
-        This returns a list of numpy or cupy arrays containing the particle array data
+        This returns a list of numpy or cupy arrays containing the particle int array data
         on each tile for this process.
 
         Unless copy_to_host is specified, the data for the arrays are not
@@ -266,12 +264,52 @@ def get_particle_arrays(self, comp_name, level, copy_to_host=False):
         List of arrays
             The requested particle array data
         '''
-        comp_idx = self.particle_container.get_comp_index(comp_name)
+        comp_idx = self.particle_container.get_icomp_index(comp_name)
 
         data_array = []
         for pti in libwarpx.libwarpx_so.WarpXParIter(self.particle_container, level):
             soa = pti.soa()
-            idx = soa.GetRealData(comp_idx)
+            idx = soa.GetIntData(comp_idx)
+            if copy_to_host:
+                data_array.append(idx.to_numpy(copy=True))
+            else:
+                xp, cupy_status = load_cupy()
+                if cupy_status is not None:
+                    libwarpx.amr.Print(cupy_status)
+
+                data_array.append(xp.array(idx, copy=False))
+
+        return data_array
+
+
+    def get_particle_idcpu_arrays(self, level, copy_to_host=False):
+        '''
+        This returns a list of numpy or cupy arrays containing the particle idcpu data
+        on each tile for this process.
+
+        Unless copy_to_host is specified, the data for the arrays are not
+        copied, but share the underlying memory buffer with WarpX. The
+        arrays are fully writeable.
+
+        Parameters
+        ----------
+        level          : int
+            The refinement level to reference (default=0)
+
+        copy_to_host   : bool
+            For GPU-enabled runs, one can either return the GPU
+            arrays (the default) or force a device-to-host copy.
+
+        Returns
+        -------
+
+        List of arrays
+            The requested particle array data
+        '''
+        data_array = []
+        for pti in libwarpx.libwarpx_so.WarpXParIter(self.particle_container, level):
+            soa = pti.soa()
+            idx = soa.GetIdCPUData()
             if copy_to_host:
                 data_array.append(idx.to_numpy(copy=True))
             else:
@@ -284,6 +322,31 @@ def get_particle_arrays(self, comp_name, level, copy_to_host=False):
         return data_array
 
 
+    def get_particle_idcpu(self, level=0, copy_to_host=False):
+        '''
+        Return a list of numpy or cupy arrays containing the particle 'idcpu'
+        numbers on each tile.
+
+        Parameters
+        ----------
+
+        level        : int
+            The refinement level to reference (default=0)
+
+        copy_to_host : bool
+            For GPU-enabled runs, one can either return the GPU
+            arrays (the default) or force a device-to-host copy.
+
+        Returns
+        -------
+
+        List of arrays
+            The requested particle idcpu
+        '''
+        return self.get_particle_idcpu_arrays(level, copy_to_host=copy_to_host)
+    idcpu = property(get_particle_idcpu)
+
+
     def get_particle_id(self, level=0, copy_to_host=False):
         '''
         Return a list of numpy or cupy arrays containing the particle 'id'
@@ -305,8 +368,8 @@ def get_particle_id(self, level=0, copy_to_host=False):
         List of arrays
             The requested particle ids
         '''
-        structs = self.get_particle_structs(level, copy_to_host)
-        return [libwarpx.amr.unpack_ids(struct['cpuid']) for struct in structs]
+        idcpu = self.get_particle_idcpu(level, copy_to_host)
+        return [libwarpx.amr.unpack_ids(tile) for tile in idcpu]
 
 
     def get_particle_cpu(self, level=0, copy_to_host=False):
@@ -330,8 +393,8 @@ def get_particle_cpu(self, level=0, copy_to_host=False):
         List of arrays
             The requested particle cpus
         '''
-        structs = self.get_particle_structs(level, copy_to_host)
-        return [libwarpx.amr.unpack_cpus(struct['cpuid']) for struct in structs]
+        idcpu = self.get_particle_idcpu(level, copy_to_host)
+        return [libwarpx.amr.unpack_cpus(tile) for tile in idcpu]
 
 
     def get_particle_x(self, level=0, copy_to_host=False):
@@ -355,20 +418,7 @@ def get_particle_x(self, level=0, copy_to_host=False):
         List of arrays
             The requested particle x position
         '''
-        if copy_to_host:
-            # Use the numpy version of cosine
-            xp = np
-        else:
-            xp, cupy_status = load_cupy()
-
-        structs = self.get_particle_structs(level, copy_to_host)
-        if libwarpx.geometry_dim == '3d' or libwarpx.geometry_dim == '2d':
-            return [struct['x'] for struct in structs]
-        elif libwarpx.geometry_dim == 'rz':
-            theta = self.get_particle_theta(level, copy_to_host)
-            return [struct['x']*xp.cos(theta) for struct, theta in zip(structs, theta)]
-        elif libwarpx.geometry_dim == '1d':
-            raise Exception('get_particle_x: There is no x coordinate with 1D Cartesian')
+        return self.get_particle_real_arrays('x', level, copy_to_host=copy_to_host)
     xp = property(get_particle_x)
 
 
@@ -393,20 +443,7 @@ def get_particle_y(self, level=0, copy_to_host=False):
         List of arrays
             The requested particle y position
         '''
-        if copy_to_host:
-            # Use the numpy version of sine
-            xp = np
-        else:
-            xp, cupy_status = load_cupy()
-
-        structs = self.get_particle_structs(level, copy_to_host)
-        if libwarpx.geometry_dim == '3d':
-            return [struct['y'] for struct in structs]
-        elif libwarpx.geometry_dim == 'rz':
-            theta = self.get_particle_theta(level, copy_to_host)
-            return [struct['x']*xp.sin(theta) for struct, theta in zip(structs, theta)]
-        elif libwarpx.geometry_dim == '1d' or libwarpx.geometry_dim == '2d':
-            raise Exception('get_particle_y: There is no y coordinate with 1D or 2D Cartesian')
+        return self.get_particle_real_arrays('y', level, copy_to_host=copy_to_host)
     yp = property(get_particle_y)
 
 
@@ -433,11 +470,12 @@ def get_particle_r(self, level=0, copy_to_host=False):
         '''
         xp, cupy_status = load_cupy()
 
-        structs = self.get_particle_structs(level, copy_to_host)
         if libwarpx.geometry_dim == 'rz':
-            return [struct['x'] for struct in structs]
+            return self.get_particle_x(level, copy_to_host)
         elif libwarpx.geometry_dim == '3d':
-            return [xp.sqrt(struct['x']**2 + struct['y']**2) for struct in structs]
+            x = self.get_particle_x(level, copy_to_host)
+            y = self.get_particle_y(level, copy_to_host)
+            return xp.sqrt(x**2 + y**2)
         elif libwarpx.geometry_dim == '2d' or libwarpx.geometry_dim == '1d':
             raise Exception('get_particle_r: There is no r coordinate with 1D or 2D Cartesian')
     rp = property(get_particle_r)
@@ -467,10 +505,11 @@ def get_particle_theta(self, level=0, copy_to_host=False):
         xp, cupy_status = load_cupy()
 
         if libwarpx.geometry_dim == 'rz':
-            return self.get_particle_arrays('theta', level, copy_to_host)
+            return self.get_particle_real_arrays('theta', level, copy_to_host)
         elif libwarpx.geometry_dim == '3d':
-            structs = self.get_particle_structs(level, copy_to_host)
-            return [xp.arctan2(struct['y'], struct['x']) for struct in structs]
+            x = self.get_particle_x(level, copy_to_host)
+            y = self.get_particle_y(level, copy_to_host)
+            return xp.arctan2(y, x)
         elif libwarpx.geometry_dim == '2d' or libwarpx.geometry_dim == '1d':
             raise Exception('get_particle_theta: There is no theta coordinate with 1D or 2D Cartesian')
     thetap = property(get_particle_theta)
@@ -497,13 +536,7 @@ def get_particle_z(self, level=0, copy_to_host=False):
         List of arrays
             The requested particle z position
         '''
-        structs = self.get_particle_structs(level, copy_to_host)
-        if libwarpx.geometry_dim == '3d':
-            return [struct['z'] for struct in structs]
-        elif libwarpx.geometry_dim == 'rz' or libwarpx.geometry_dim == '2d':
-            return [struct['y'] for struct in structs]
-        elif libwarpx.geometry_dim == '1d':
-            return [struct['x'] for struct in structs]
+        return self.get_particle_real_arrays('z', level, copy_to_host=copy_to_host)
     zp = property(get_particle_z)
 
 
@@ -528,7 +561,7 @@ def get_particle_weight(self, level=0, copy_to_host=False):
         List of arrays
             The requested particle weight
         '''
-        return self.get_particle_arrays('w', level, copy_to_host=copy_to_host)
+        return self.get_particle_real_arrays('w', level, copy_to_host=copy_to_host)
     wp = property(get_particle_weight)
 
 
@@ -553,7 +586,7 @@ def get_particle_ux(self, level=0, copy_to_host=False):
         List of arrays
             The requested particle x momentum
         '''
-        return self.get_particle_arrays('ux', level, copy_to_host=copy_to_host)
+        return self.get_particle_real_arrays('ux', level, copy_to_host=copy_to_host)
     uxp = property(get_particle_ux)
 
 
@@ -578,7 +611,7 @@ def get_particle_uy(self, level=0, copy_to_host=False):
         List of arrays
             The requested particle y momentum
         '''
-        return self.get_particle_arrays('uy', level, copy_to_host=copy_to_host)
+        return self.get_particle_real_arrays('uy', level, copy_to_host=copy_to_host)
     uyp = property(get_particle_uy)
 
 
@@ -604,7 +637,7 @@ def get_particle_uz(self, level=0, copy_to_host=False):
             The requested particle z momentum
         '''
 
-        return self.get_particle_arrays('uz', level, copy_to_host=copy_to_host)
+        return self.get_particle_real_arrays('uz', level, copy_to_host=copy_to_host)
     uzp = property(get_particle_uz)
 
 
@@ -720,70 +753,6 @@ def get_particle_boundary_buffer_size(self, species_name, boundary, local=False)
         )
 
 
-    def get_particle_boundary_buffer_structs(
-            self, species_name, boundary, level, copy_to_host=False
-        ):
-        '''
-        This returns a list of numpy or cupy arrays containing the particle struct data
-        for a species that has been scraped by a specific simulation boundary. The
-        particle data is represented as a structured array and contains the
-        particle 'x', 'y', 'z', and 'idcpu'.
-
-        Unless copy_to_host is specified, the data for the structs are
-        not copied, but share the underlying memory buffer with WarpX. The
-        arrays are fully writeable.
-
-        Note that cupy does not support structs:
-        https://github.com/cupy/cupy/issues/2031
-        and will return arrays of binary blobs for the AoS (DP: ``"|V24"``). If copied
-        to host via copy_to_host, we correct for the right numpy AoS type.
-
-        Parameters
-        ----------
-
-        species_name : str
-            The species name that the data will be returned for
-
-        boundary     : str
-            The boundary from which to get the scraped particle data in the
-            form x/y/z_hi/lo or eb.
-
-        level        : int
-            The refinement level to reference (default=0)
-
-        copy_to_host : bool
-            For GPU-enabled runs, one can either return the GPU
-            arrays (the default) or force a device-to-host copy.
-
-        Returns
-        -------
-
-        List of arrays
-            The requested particle struct data
-        '''
-        particle_container = self.particle_buffer.get_particle_container(
-            species_name, self._get_boundary_number(boundary)
-        )
-
-        particle_data = []
-        for pti in libwarpx.libwarpx_so.BoundaryBufferParIter(particle_container, level):
-            if copy_to_host:
-                particle_data.append(pti.aos().to_numpy(copy=True))
-            else:
-                if libwarpx.amr.Config.have_gpu:
-                    libwarpx.amr.Print(
-                        "get_particle_structs: cupy does not yet support structs. "
-                        "https://github.com/cupy/cupy/issues/2031"
-                        "Did you mean copy_to_host=True?"
-                    )
-                xp, cupy_status = load_cupy()
-                if cupy_status is not None:
-                    libwarpx.amr.Print(cupy_status)
-                aos_arr = xp.array(pti.aos(), copy=False)   # void blobs for cupy
-                particle_data.append(aos_arr)
-        return particle_data
-
-
     def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level):
         '''
         This returns a list of numpy or cupy arrays containing the particle array data
diff --git a/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H
index e88b522f148..fa00f6288f9 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H
+++ b/Source/Diagnostics/ComputeDiagFunctors/BackTransformParticleFunctor.H
@@ -143,19 +143,19 @@ struct LorentzTransformParticles
         const amrex::ParticleReal uzp = uz_old_p * weight_old
                                       + uz_new_p * weight_new;
 #if defined (WARPX_DIM_3D)
-        dst.m_aos[i_dst].pos(0) = xp;
-        dst.m_aos[i_dst].pos(1) = yp;
-        dst.m_aos[i_dst].pos(2) = zp;
+        dst.m_rdata[PIdx::x][i_dst] = xp;
+        dst.m_rdata[PIdx::y][i_dst] = yp;
+        dst.m_rdata[PIdx::z][i_dst] = zp;
 #elif defined (WARPX_DIM_RZ)
-        dst.m_aos[i_dst].pos(0) = std::sqrt(xp*xp + yp*yp);
-        dst.m_aos[i_dst].pos(1) = zp;
+        dst.m_rdata[PIdx::x][i_dst] = std::sqrt(xp*xp + yp*yp);
+        dst.m_rdata[PIdx::z][i_dst] = zp;
         dst.m_rdata[PIdx::theta][i_dst] = std::atan2(yp, xp);
 #elif defined (WARPX_DIM_XZ)
-        dst.m_aos[i_dst].pos(0) = xp;
-        dst.m_aos[i_dst].pos(1) = zp;
+        dst.m_rdata[PIdx::x][i_dst] = xp;
+        dst.m_rdata[PIdx::z][i_dst] = zp;
         amrex::ignore_unused(yp);
 #elif defined (WARPX_DIM_1D_Z)
-        dst.m_aos[i_dst].pos(0) = zp;
+        dst.m_rdata[PIdx::z][i_dst] = zp;
         amrex::ignore_unused(xp, yp);
 #else
         amrex::ignore_unused(xp, yp, zp);
diff --git a/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp b/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp
index abfba37cd15..0024122ebf9 100644
--- a/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp
+++ b/Source/Diagnostics/FlushFormats/FlushFormatAscent.cpp
@@ -94,9 +94,6 @@ FlushFormatAscent::WriteParticles(const amrex::Vector<ParticleDiag>& particle_di
         // get names of real comps
         std::map<std::string, int> real_comps_map = pc->getParticleComps();
 
-        // WarpXParticleContainer compile-time extra AoS attributes (Real): 0
-        // WarpXParticleContainer compile-time extra AoS attributes (int): 0
-
         // WarpXParticleContainer compile-time extra SoA attributes (Real): PIdx::nattribs
         // not an efficient search, but N is small...
         for(int j = 0; j < PIdx::nattribs; ++j)
diff --git a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp
index d77437fb931..b083e60529f 100644
--- a/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp
+++ b/Source/Diagnostics/FlushFormats/FlushFormatCheckpoint.cpp
@@ -178,8 +178,8 @@ FlushFormatCheckpoint::CheckpointParticles (
         Vector<std::string> real_names;
         Vector<std::string> int_names;
 
+        // note: positions skipped here, since we reconstruct a plotfile SoA from them
         real_names.push_back("weight");
-
         real_names.push_back("momentum_x");
         real_names.push_back("momentum_y");
         real_names.push_back("momentum_z");
@@ -189,9 +189,12 @@ FlushFormatCheckpoint::CheckpointParticles (
 #endif
 
         // get the names of the real comps
-        real_names.resize(pc->NumRealComps());
+        //   note: skips the mandatory AMREX_SPACEDIM positions for pure SoA
+        real_names.resize(pc->NumRealComps() - AMREX_SPACEDIM);
         auto runtime_rnames = pc->getParticleRuntimeComps();
-        for (auto const& x : runtime_rnames) { real_names[x.second+PIdx::nattribs] = x.first; }
+        for (auto const& x : runtime_rnames) {
+            real_names[x.second + PIdx::nattribs - AMREX_SPACEDIM] = x.first;
+        }
 
         // and the int comps
         int_names.resize(pc->NumIntComps());
diff --git a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp
index 970d9a504d2..880e2df01ff 100644
--- a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp
+++ b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp
@@ -355,8 +355,8 @@ FlushFormatPlotfile::WriteParticles(const std::string& dir,
         Vector<int> int_flags;
         Vector<int> real_flags;
 
+        // note: positions skipped here, since we reconstruct a plotfile SoA from them
         real_names.push_back("weight");
-
         real_names.push_back("momentum_x");
         real_names.push_back("momentum_y");
         real_names.push_back("momentum_z");
@@ -366,14 +366,21 @@ FlushFormatPlotfile::WriteParticles(const std::string& dir,
 #endif
 
         // get the names of the real comps
-        real_names.resize(tmp.NumRealComps());
+
+        //   note: skips the mandatory AMREX_SPACEDIM positions for pure SoA
+        real_names.resize(tmp.NumRealComps() - AMREX_SPACEDIM);
         auto runtime_rnames = tmp.getParticleRuntimeComps();
-        for (auto const& x : runtime_rnames) { real_names[x.second+PIdx::nattribs] = x.first; }
+        for (auto const& x : runtime_rnames) {
+            real_names[x.second + PIdx::nattribs - AMREX_SPACEDIM] = x.first;
+        }
 
         // plot any "extra" fields by default
         real_flags = part_diag.m_plot_flags;
         real_flags.resize(tmp.NumRealComps(), 1);
 
+        //   note: skip the mandatory AMREX_SPACEDIM positions for pure SoA
+        real_flags.erase(real_flags.begin(), real_flags.begin() + AMREX_SPACEDIM);
+
         // and the names
         int_names.resize(tmp.NumIntComps());
         auto runtime_inames = tmp.getParticleRuntimeiComps();
diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp
index 7ca5e6541d7..a8bb9303fe1 100644
--- a/Source/Diagnostics/ParticleIO.cpp
+++ b/Source/Diagnostics/ParticleIO.cpp
@@ -160,7 +160,7 @@ MultiParticleContainer::Restart (const std::string& dir)
             );
         }
 
-        for (int j = PIdx::nattribs; j < nr; ++j) {
+        for (int j = PIdx::nattribs-AMREX_SPACEDIM; j < nr; ++j) {
             const auto& comp_name = real_comp_names[j];
             auto current_comp_names = pc->getParticleComps();
             auto search = current_comp_names.find(comp_name);
diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp
index 9f45392bb0a..24ad0e64ea8 100644
--- a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp
+++ b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp
@@ -431,8 +431,6 @@ void FieldProbe::ComputeDiags (int step)
         {
             const auto getPosition = GetParticlePosition<FieldProbePIdx>(pti);
             auto setPosition = SetParticlePosition<FieldProbePIdx>(pti);
-            const auto& aos = pti.GetArrayOfStructs();
-            const auto* AMREX_RESTRICT m_structs = aos().dataPtr();
 
             auto const np = pti.numParticles();
             if (update_particles_moving_window)
@@ -482,6 +480,8 @@ void FieldProbe::ComputeDiags (int step)
                 ParticleReal* const AMREX_RESTRICT part_Bz = attribs[FieldProbePIdx::Bz].dataPtr();
                 ParticleReal* const AMREX_RESTRICT part_S = attribs[FieldProbePIdx::S].dataPtr();
 
+                auto * const AMREX_RESTRICT idcpu = pti.GetStructOfArrays().GetIdCPUData().data();
+
                 const auto &xyzmin = WarpX::LowerCorner(box, lev, 0._rt);
                 const std::array<Real, 3> &dx = WarpX::CellSize(lev);
 
@@ -556,7 +556,7 @@ void FieldProbe::ComputeDiags (int step)
                         amrex::ParticleReal xp, yp, zp;
                         getPosition(ip, xp, yp, zp);
                         long idx = ip*noutputs;
-                        dvp[idx++] = m_structs[ip].id();
+                        dvp[idx++] = amrex::ParticleIDWrapper{idcpu[ip]};  // all particles created on IO cpu
                         dvp[idx++] = xp;
                         dvp[idx++] = yp;
                         dvp[idx++] = zp;
diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H
index c85bf8fd541..7d59ade5dc6 100644
--- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H
+++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H
@@ -24,7 +24,14 @@ struct FieldProbePIdx
 {
     enum
     {
-        Ex = 0, Ey, Ez,
+#if !defined (WARPX_DIM_1D_Z)
+        x,
+#endif
+#if defined (WARPX_DIM_3D)
+        y,
+#endif
+        z,
+        Ex, Ey, Ez,
         Bx, By, Bz,
         S, //!< the Poynting vector
 #ifdef WARPX_DIM_RZ
@@ -40,9 +47,14 @@ struct FieldProbePIdx
  * nattribs tells the particle container to allot 7 SOA values.
  */
 class FieldProbeParticleContainer
-    : public amrex::ParticleContainer<0, 0, FieldProbePIdx::nattribs>
+    : public amrex::ParticleContainerPureSoA<FieldProbePIdx::nattribs, 0>
 {
 public:
+    static constexpr int NStructReal = 0;
+    static constexpr int NStructInt = 0;
+    static constexpr int NReal = FieldProbePIdx::nattribs;
+    static constexpr int NInt = 0;
+
     FieldProbeParticleContainer (amrex::AmrCore* amr_core);
     ~FieldProbeParticleContainer() override = default;
 
@@ -52,9 +64,9 @@ public:
     FieldProbeParticleContainer& operator= ( FieldProbeParticleContainer&& )       = default;
 
     //! amrex iterator for our number of attributes
-    using iterator = amrex::ParIter<0, 0, FieldProbePIdx::nattribs, 0>;
+    using iterator = amrex::ParIterSoA<FieldProbePIdx::nattribs, 0>;
     //! amrex iterator for our number of attributes (read-only)
-    using const_iterator = amrex::ParConstIter<0, 0, FieldProbePIdx::nattribs, 0>;
+    using const_iterator = amrex::ParConstIterSoA<FieldProbePIdx::nattribs, 0>;
 
     //! similar to WarpXParticleContainer::AddNParticles but does not include u(x,y,z)
     void AddNParticles (int lev, amrex::Vector<amrex::ParticleReal> const & x, amrex::Vector<amrex::ParticleReal> const & y, amrex::Vector<amrex::ParticleReal> const & z);
diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp
index 1fd741ddc47..18137fe9b2c 100644
--- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp
+++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp
@@ -59,7 +59,7 @@
 using namespace amrex;
 
 FieldProbeParticleContainer::FieldProbeParticleContainer (AmrCore* amr_core)
-    : ParticleContainer<0, 0, FieldProbePIdx::nattribs>(amr_core->GetParGDB())
+    : ParticleContainerPureSoA<FieldProbePIdx::nattribs, 0>(amr_core->GetParGDB())
 {
     SetParticleSize();
 }
@@ -89,33 +89,17 @@ FieldProbeParticleContainer::AddNParticles (int lev,
      * is then coppied to the permament tile which is stored on the particle
      * (particle_tile).
      */
+    using PinnedTile = typename ContainerLike<amrex::PinnedArenaAllocator>::ParticleTileType;
 
-    using PinnedTile = ParticleTile<amrex::Particle<NStructReal, NStructInt>,
-                                    NArrayReal, NArrayInt,
-                                    amrex::PinnedArenaAllocator>;
     PinnedTile pinned_tile;
     pinned_tile.define(NumRuntimeRealComps(), NumRuntimeIntComps());
 
     for (int i = 0; i < np; i++)
     {
-        ParticleType p;
-        p.id() = ParticleType::NextID();
-        p.cpu() = ParallelDescriptor::MyProc();
-#if defined(WARPX_DIM_3D)
-        p.pos(0) = x[i];
-        p.pos(1) = y[i];
-        p.pos(2) = z[i];
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-        amrex::ignore_unused(y);
-        p.pos(0) = x[i];
-        p.pos(1) = z[i];
-#elif defined(WARPX_DIM_1D_Z)
-        amrex::ignore_unused(x, y);
-        p.pos(0) = z[i];
-#endif
-
-        // write position, cpu id, and particle id to particle
-        pinned_tile.push_back(p);
+        auto & idcpu_data = pinned_tile.GetStructOfArrays().GetIdCPUData();
+        idcpu_data.push_back(0);
+        amrex::ParticleIDWrapper{idcpu_data.back()} = ParticleType::NextID();
+        amrex::ParticleCPUWrapper(idcpu_data.back()) = ParallelDescriptor::MyProc();
     }
 
     // write Real attributes (SoA) to particle initialized zero
@@ -125,7 +109,13 @@ FieldProbeParticleContainer::AddNParticles (int lev,
 #ifdef WARPX_DIM_RZ
     pinned_tile.push_back_real(FieldProbePIdx::theta, np, 0.0);
 #endif
-
+#if !defined (WARPX_DIM_1D_Z)
+    pinned_tile.push_back_real(FieldProbePIdx::x, x);
+#endif
+#if defined (WARPX_DIM_3D)
+    pinned_tile.push_back_real(FieldProbePIdx::y, y);
+#endif
+    pinned_tile.push_back_real(FieldProbePIdx::z, z);
     pinned_tile.push_back_real(FieldProbePIdx::Ex, np, 0.0);
     pinned_tile.push_back_real(FieldProbePIdx::Ey, np, 0.0);
     pinned_tile.push_back_real(FieldProbePIdx::Ez, np, 0.0);
diff --git a/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp b/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp
index 893b00a5f00..b4e07b51982 100644
--- a/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp
+++ b/Source/Diagnostics/ReducedDiags/LoadBalanceCosts.cpp
@@ -56,8 +56,7 @@ namespace
             auto const & plev  = pc.GetParticles(lev);
 
             auto const & ptile = plev.at(box_index);
-            auto const & aos   = ptile.GetArrayOfStructs();
-            auto const np = aos.numParticles();
+            auto const np = ptile.numParticles();
             num_macro_particles += np;
         }
 
diff --git a/Source/Diagnostics/WarpXOpenPMD.H b/Source/Diagnostics/WarpXOpenPMD.H
index 4597dacd9ae..6c904790e15 100644
--- a/Source/Diagnostics/WarpXOpenPMD.H
+++ b/Source/Diagnostics/WarpXOpenPMD.H
@@ -41,7 +41,7 @@ class WarpXParticleCounter
 {
 public:
   using ParticleContainer = typename WarpXParticleContainer::ContainerLike<amrex::PinnedArenaAllocator>;
-  using ParticleIter = typename amrex::ParIter<0, 0, PIdx::nattribs, 0, amrex::PinnedArenaAllocator>;
+  using ParticleIter = typename amrex::ParIterSoA<PIdx::nattribs, 0, amrex::PinnedArenaAllocator>;
 
   WarpXParticleCounter (ParticleContainer* pc);
   [[nodiscard]] unsigned long GetTotalNumParticles () const {return m_Total;}
@@ -77,7 +77,7 @@ class WarpXOpenPMDPlot
 {
 public:
   using ParticleContainer = typename WarpXParticleContainer::ContainerLike<amrex::PinnedArenaAllocator>;
-  using ParticleIter = typename amrex::ParConstIter<0, 0, PIdx::nattribs, 0, amrex::PinnedArenaAllocator>;
+  using ParticleIter = typename amrex::ParConstIterSoA<PIdx::nattribs, 0, amrex::PinnedArenaAllocator>;
 
   /** Initialize openPMD I/O routines
    *
diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp
index 7cc9f571a4a..39717ef6ec5 100644
--- a/Source/Diagnostics/WarpXOpenPMD.cpp
+++ b/Source/Diagnostics/WarpXOpenPMD.cpp
@@ -18,11 +18,9 @@
 #include "WarpX.H"
 #include "OpenPMDHelpFunction.H"
 
-#include <ablastr/particles/IndexHandling.H>
 #include <ablastr/warn_manager/WarnManager.H>
 
 #include <AMReX.H>
-#include <AMReX_ArrayOfStructs.H>
 #include <AMReX_BLassert.H>
 #include <AMReX_Box.H>
 #include <AMReX_Config.H>
@@ -550,6 +548,13 @@ for (unsigned i = 0, n = particle_diags.size(); i < n; ++i) {
     // see openPMD ED-PIC extension for namings
     // note: an underscore separates the record name from its component
     //       for non-scalar records
+#if !defined (WARPX_DIM_1D_Z)
+    real_names.push_back("position_x");
+#endif
+#if defined (WARPX_DIM_3D)
+    real_names.push_back("position_y");
+#endif
+    real_names.push_back("position_z");
     real_names.push_back("weighting");
     real_names.push_back("momentum_x");
     real_names.push_back("momentum_y");
@@ -722,77 +727,7 @@ WarpXOpenPMDPlot::DumpToFile (ParticleContainer* pc,
 
             contributed_particles = true;
 
-            // get position and particle ID from aos
-            // note: this implementation iterates the AoS 4x...
-            // if we flush late as we do now, we can also copy out the data in one go
-            const auto &aos = pti.GetArrayOfStructs();  // size =  numParticlesOnTile
-            {
-                // Save positions
-#if defined(WARPX_DIM_RZ)
-                {
-                   const std::shared_ptr<amrex::ParticleReal> z(
-                           new amrex::ParticleReal[numParticleOnTile],
-                           [](amrex::ParticleReal const *p) { delete[] p; }
-                   );
-                   for (auto i = 0; i < numParticleOnTile; i++) {
-                       z.get()[i] = aos[i].pos(1);  // {0: "r", 1: "z"}
-                   }
-                   std::string const positionComponent = "z";
-                   currSpecies["position"]["z"].storeChunk(z, {offset}, {numParticleOnTile64});
-                }
-
-                //   reconstruct x and y from polar coordinates r, theta
-                auto const& soa = pti.GetStructOfArrays();
-                amrex::ParticleReal const* theta = soa.GetRealData(PIdx::theta).dataPtr();
-                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(theta != nullptr, "openPMD: invalid theta pointer.");
-                WARPX_ALWAYS_ASSERT_WITH_MESSAGE(int(soa.GetRealData(PIdx::theta).size()) == numParticleOnTile,
-                                                 "openPMD: theta and tile size do not match");
-                {
-                    const std::shared_ptr< amrex::ParticleReal > x(
-                            new amrex::ParticleReal[numParticleOnTile],
-                            [](amrex::ParticleReal const *p){ delete[] p; }
-                    );
-                    const std::shared_ptr< amrex::ParticleReal > y(
-                            new amrex::ParticleReal[numParticleOnTile],
-                            [](amrex::ParticleReal const *p){ delete[] p; }
-                    );
-                    for (auto i=0; i<numParticleOnTile; i++) {
-                        auto const r = aos[i].pos(0);  // {0: "r", 1: "z"}
-                        x.get()[i] = r * std::cos(theta[i]);
-                        y.get()[i] = r * std::sin(theta[i]);
-                    }
-                    currSpecies["position"]["x"].storeChunk(x, {offset}, {numParticleOnTile64});
-                    currSpecies["position"]["y"].storeChunk(y, {offset}, {numParticleOnTile64});
-                }
-#else
-                auto const positionComponents = detail::getParticlePositionComponentLabels();
-                for (auto currDim = 0; currDim < AMREX_SPACEDIM; currDim++) {
-                    const std::shared_ptr<amrex::ParticleReal> curr(
-                            new amrex::ParticleReal[numParticleOnTile],
-                            [](amrex::ParticleReal const *p) { delete[] p; }
-                    );
-                    for (auto i = 0; i < numParticleOnTile; i++) {
-                        curr.get()[i] = aos[i].pos(currDim);
-                    }
-                    std::string const positionComponent = positionComponents[currDim];
-                    currSpecies["position"][positionComponent].storeChunk(curr, {offset},
-                                                                          {numParticleOnTile64});
-                }
-#endif
-
-                // save particle ID after converting it to a globally unique ID
-                const std::shared_ptr<uint64_t> ids(
-                        new uint64_t[numParticleOnTile],
-                        [](uint64_t const *p) { delete[] p; }
-                );
-                for (auto i = 0; i < numParticleOnTile; i++) {
-                    ids.get()[i] = ablastr::particles::localIDtoGlobal(static_cast<int>(aos[i].id()), static_cast<int>(aos[i].cpu()));
-                }
-                const auto *const scalar = openPMD::RecordComponent::SCALAR;
-                currSpecies["id"][scalar].storeChunk(ids, {offset}, {numParticleOnTile64});
-
-            }
-            //  save "extra" particle properties in AoS and SoA
+            //  save particle properties
             SaveRealProperty(pti,
                              currSpecies,
                              offset,
@@ -893,10 +828,9 @@ WarpXOpenPMDPlot::SetupRealProperties (ParticleContainer const * pc,
 
     std::set< std::string > addedRecords; // add meta-data per record only once
     for (auto idx=0; idx<pc->NumRealComps(); idx++) {
-        auto ii = ParticleContainer::NStructReal + idx; // jump over extra AoS names
-        if (write_real_comp[ii]) {
+        if (write_real_comp[idx]) {
             // handle scalar and non-scalar records by name
-            const auto [record_name, component_name] = detail::name2openPMD(real_comp_names[ii]);
+            const auto [record_name, component_name] = detail::name2openPMD(real_comp_names[idx]);
             auto currRecord = currSpecies[record_name];
 
             // meta data for ED-PIC extension
@@ -917,10 +851,9 @@ WarpXOpenPMDPlot::SetupRealProperties (ParticleContainer const * pc,
         }
     }
     for (auto idx=0; idx<int_counter; idx++) {
-        auto ii = ParticleContainer::NStructInt + idx; // jump over extra AoS names
-        if (write_int_comp[ii]) {
+        if (write_int_comp[idx]) {
             // handle scalar and non-scalar records by name
-            const auto [record_name, component_name] = detail::name2openPMD(int_comp_names[ii]);
+            const auto [record_name, component_name] = detail::name2openPMD(int_comp_names[idx]);
             auto currRecord = currSpecies[record_name];
 
             // meta data for ED-PIC extension
@@ -949,33 +882,8 @@ WarpXOpenPMDPlot::SaveRealProperty (ParticleIter& pti,
 
 {
   auto const numParticleOnTile = pti.numParticles();
-  auto const numParticleOnTile64 = static_cast<uint64_t>( numParticleOnTile );
-  auto const& aos = pti.GetArrayOfStructs();  // size =  numParticlesOnTile
+  auto const numParticleOnTile64 = static_cast<uint64_t>(numParticleOnTile);
   auto const& soa = pti.GetStructOfArrays();
-  // first we concatenate the AoS into contiguous arrays
-  {
-    // note: WarpX does not yet use extra AoS Real attributes
-    for( auto idx=0; idx<ParticleIter::ContainerType::NStructReal; idx++ ) {  // lgtm [cpp/constant-comparison]
-        if( write_real_comp[idx] ) {
-            // handle scalar and non-scalar records by name
-            const auto [record_name, component_name] = detail::name2openPMD(real_comp_names[idx]);
-            auto currRecord = currSpecies[record_name];
-            auto currRecordComp = currRecord[component_name];
-
-            const std::shared_ptr< amrex::ParticleReal > d(
-                new amrex::ParticleReal[numParticleOnTile],
-                [](amrex::ParticleReal const *p){ delete[] p; }
-            );
-
-            for( auto kk=0; kk<numParticleOnTile; kk++ ) {
-                d.get()[kk] = aos[kk].rdata(idx);
-            }
-
-            currRecordComp.storeChunk(d,
-                {offset}, {numParticleOnTile64});
-        }
-    }
-  }
 
   auto const getComponentRecord = [&currSpecies](std::string const comp_name) {
     // handle scalar and non-scalar records by name
@@ -983,24 +891,29 @@ WarpXOpenPMDPlot::SaveRealProperty (ParticleIter& pti,
     return currSpecies[record_name][component_name];
   };
 
+    // here we the save the SoA properties (idcpu)
+    {
+        // todo: add support to not write the particle index
+        getComponentRecord("id").storeChunkRaw(
+        soa.GetIdCPUData().data(), {offset}, {numParticleOnTile64});
+    }
+
   // here we the save the SoA properties (real)
   {
     auto const real_counter = std::min(write_real_comp.size(), real_comp_names.size());
     for (auto idx=0; idx<real_counter; idx++) {
-        auto ii = ParticleIter::ContainerType::NStructReal + idx;  // jump over extra AoS names
-        if (write_real_comp[ii]) {
-            getComponentRecord(real_comp_names[ii]).storeChunkRaw(
-                soa.GetRealData(idx).data(), {offset}, {numParticleOnTile64});
-        }
+      if (write_real_comp[idx]) {
+        getComponentRecord(real_comp_names[idx]).storeChunkRaw(
+          soa.GetRealData(idx).data(), {offset}, {numParticleOnTile64});
+      }
     }
   }
   // and now SoA int properties
   {
     auto const int_counter = std::min(write_int_comp.size(), int_comp_names.size());
     for (auto idx=0; idx<int_counter; idx++) {
-        auto ii = ParticleIter::ContainerType::NStructInt + idx;  // jump over extra AoS names
-        if (write_int_comp[ii]) {
-            getComponentRecord(int_comp_names[ii]).storeChunkRaw(
+        if (write_int_comp[idx]) {
+            getComponentRecord(int_comp_names[idx]).storeChunkRaw(
                 soa.GetIntData(idx).data(), {offset}, {numParticleOnTile64});
         }
     }
diff --git a/Source/EmbeddedBoundary/ParticleBoundaryProcess.H b/Source/EmbeddedBoundary/ParticleBoundaryProcess.H
index a78aa19458a..2c8211fadf6 100644
--- a/Source/EmbeddedBoundary/ParticleBoundaryProcess.H
+++ b/Source/EmbeddedBoundary/ParticleBoundaryProcess.H
@@ -25,12 +25,11 @@ struct NoOp {
 struct Absorb {
     template <typename PData>
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    void operator() (const PData& ptd, int i,
+    void operator() (PData& ptd, int i,
                      const amrex::RealVect& /*pos*/, const amrex::RealVect& /*normal*/,
                      amrex::RandomEngine const& /*engine*/) const noexcept
     {
-        auto& p = ptd.m_aos[i];
-        p.id() = -p.id();
+        ptd.id(i) = -ptd.id(i);
     }
 };
 }
diff --git a/Source/EmbeddedBoundary/ParticleScraper.H b/Source/EmbeddedBoundary/ParticleScraper.H
index d6196c35f44..312b3762a7e 100644
--- a/Source/EmbeddedBoundary/ParticleScraper.H
+++ b/Source/EmbeddedBoundary/ParticleScraper.H
@@ -170,13 +170,12 @@ scrapeParticles (PC& pc, const amrex::Vector<const amrex::MultiFab*>& distance_t
             auto& tile = pti.GetParticleTile();
             auto ptd = tile.getParticleTileData();
             const auto np = tile.numParticles();
-            amrex::Particle<0,0> * const particles = tile.GetArrayOfStructs()().data();
             auto phi = (*distance_to_eb[lev])[pti].array();  // signed distance function
             amrex::ParallelForRNG( np,
             [=] AMREX_GPU_DEVICE (const int ip, amrex::RandomEngine const& engine) noexcept
             {
                 // skip particles that are already flagged for removal
-                if (particles[ip].id() < 0) return;
+                if (ptd.id(ip) < 0) return;
 
                 amrex::ParticleReal xp, yp, zp;
                 getPosition(ip, xp, yp, zp);
diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
index 5c90dab25e6..c69f07acdb2 100644
--- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
+++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H
@@ -72,7 +72,8 @@ class BinaryCollision final
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
     using ParticleTileType = WarpXParticleContainer::ParticleTileType;
-    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using ParticleTileDataType = ParticleTileType::ParticleTileDataType;
+    using ParticleBins = amrex::DenseBins<ParticleTileDataType>;
     using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
     using index_type = ParticleBins::index_type;
 
@@ -261,9 +262,6 @@ public:
             const amrex::ParticleReal q1 = species_1.getCharge();
             const amrex::ParticleReal m1 = species_1.getMass();
             auto get_position_1  = GetParticlePosition<PIdx>(ptile_1, getpos_offset);
-            // Needed to access the particle id
-            ParticleType * AMREX_RESTRICT
-                                 particle_ptr_1 = ptile_1.GetArrayOfStructs()().data();
 
             amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev);
 #if defined WARPX_DIM_1D_Z
@@ -371,7 +369,7 @@ public:
                                                     soa_1, soa_1,
                                                     product_species_vector,
                                                     tile_products_data,
-                                                    particle_ptr_1, particle_ptr_1, m1, m1,
+                                                    m1, m1,
                                                     products_mass, p_mask, products_np,
                                                     copy_species1, copy_species2,
                                                     p_pair_indices_1, p_pair_indices_2,
@@ -403,9 +401,6 @@ public:
             const amrex::ParticleReal q1 = species_1.getCharge();
             const amrex::ParticleReal m1 = species_1.getMass();
             auto get_position_1  = GetParticlePosition<PIdx>(ptile_1, getpos_offset);
-            // Needed to access the particle id
-            ParticleType * AMREX_RESTRICT
-                                 particle_ptr_1 = ptile_1.GetArrayOfStructs()().data();
             // - Species 2
             const auto soa_2 = ptile_2.getParticleTileData();
             index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr();
@@ -413,9 +408,6 @@ public:
             const amrex::ParticleReal q2 = species_2.getCharge();
             const amrex::ParticleReal m2 = species_2.getMass();
             auto get_position_2  = GetParticlePosition<PIdx>(ptile_2, getpos_offset);
-            // Needed to access the particle id
-            ParticleType * AMREX_RESTRICT
-                                 particle_ptr_2 = ptile_2.GetArrayOfStructs()().data();
 
             amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev);
 #if defined WARPX_DIM_1D_Z
@@ -535,7 +527,7 @@ public:
                                                     soa_1, soa_2,
                                                     product_species_vector,
                                                     tile_products_data,
-                                                    particle_ptr_1, particle_ptr_2, m1, m2,
+                                                    m1, m2,
                                                     products_mass, p_mask, products_np,
                                                     copy_species1, copy_species2,
                                                     p_pair_indices_1, p_pair_indices_2,
diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H
index feb7acf81d3..cfdc36d3c50 100644
--- a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H
@@ -23,10 +23,13 @@
  * \brief This functor performs pairwise Coulomb collision on a single cell by calling the function
  *        ElasticCollisionPerez. It also reads and contains the Coulomb logarithm.
  */
-class PairWiseCoulombCollisionFunc{
+class PairWiseCoulombCollisionFunc
+{
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
-    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using ParticleTileType = WarpXParticleContainer::ParticleTileType;
+    using ParticleTileDataType = ParticleTileType::ParticleTileDataType;
+    using ParticleBins = amrex::DenseBins<ParticleTileDataType>;
     using index_type = ParticleBins::index_type;
     using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
 
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H
index c1be307b811..ab01eba2c81 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMC.H
@@ -38,7 +38,8 @@ class DSMC final
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
     using ParticleTileType = WarpXParticleContainer::ParticleTileType;
-    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using ParticleTileDataType = ParticleTileType::ParticleTileDataType;
+    using ParticleBins = amrex::DenseBins<ParticleTileDataType>;
     using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
     using index_type = ParticleBins::index_type;
 
diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
index c1fb7ee7e38..5a7a4f3237a 100644
--- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H
@@ -10,6 +10,9 @@
 #define SPLIT_AND_SCATTER_FUNC_H_
 
 #include "Particles/Collision/ScatteringProcess.H"
+#include "Particles/NamedComponentParticleContainer.H"
+
+#include <AMReX_Particle.H>
 
 /**
  * \brief Function that performs the particle scattering and injection due
@@ -55,8 +58,6 @@ int splitScatteringParticles (
     const auto ptile1_data = ptile1.getParticleTileData();
     const auto ptile2_data = ptile2.getParticleTileData();
 
-    const Long minus_one_long = -1;
-
     ParallelForRNG(n_total_pairs,
     [=] AMREX_GPU_DEVICE (int i, RandomEngine const& engine) noexcept
     {
@@ -70,20 +71,37 @@ int splitScatteringParticles (
             // starting with the parent particles
             auto& w1 = ptile1_data.m_rdata[PIdx::w][p_pair_indices_1[i]];
             auto& w2 = ptile2_data.m_rdata[PIdx::w][p_pair_indices_2[i]];
+            uint64_t* AMREX_RESTRICT idcpu1 = ptile1_data.m_idcpu;
+            uint64_t* AMREX_RESTRICT idcpu2 = ptile2_data.m_idcpu;
+
+            // Note: Particle::atomicSetID should also be provided as a standalone helper function in AMReX
+            //       to replace the following lambda.
+            auto const atomicSetIdMinus = [] AMREX_GPU_DEVICE (uint64_t & idcpu)
+            {
+                constexpr amrex::Long minus_one_long = -1;
+                uint64_t tmp = 0;
+                amrex::ParticleIDWrapper wrapper(tmp);
+                wrapper = minus_one_long;
+#if defined(AMREX_USE_OMP)
+#pragma omp atomic write
+                idcpu = wrapper.m_idata;
+#else
+                auto *old_ptr = reinterpret_cast<unsigned long long*>(&idcpu);
+                    amrex::Gpu::Atomic::Exch(old_ptr, (unsigned long long) wrapper.m_idata);
+#endif
+            };
 
             // Remove p_pair_reaction_weight[i] from the colliding particles' weights.
             // If the colliding particle weight decreases to zero, remove particle by
             // setting its id to -1.
             Gpu::Atomic::AddNoRet(&w1, -p_pair_reaction_weight[i]);
             if (w1 <= 0._prt) {
-                auto& p = ptile1_data.m_aos[p_pair_indices_1[i]];
-                p.atomicSetID(minus_one_long);
+                atomicSetIdMinus(idcpu1[p_pair_indices_1[i]]);
             }
 
             Gpu::Atomic::AddNoRet(&w2, -p_pair_reaction_weight[i]);
             if (w2 <= 0._prt) {
-                auto& p = ptile2_data.m_aos[p_pair_indices_2[i]];
-                p.atomicSetID(minus_one_long);
+                atomicSetIdMinus(idcpu2[p_pair_indices_2[i]]);
             }
 
             // Set the child particle properties appropriately
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
index 397536b67bf..b2a2112ca68 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H
@@ -33,10 +33,13 @@
  *  creation functor.
  *  This functor also reads and contains the fusion multiplier.
  */
-class NuclearFusionFunc{
+class NuclearFusionFunc
+{
     // Define shortcuts for frequently-used type names
     using ParticleType = WarpXParticleContainer::ParticleType;
-    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using ParticleTileType = WarpXParticleContainer::ParticleTileType;
+    using ParticleTileDataType = ParticleTileType::ParticleTileDataType;
+    using ParticleBins = amrex::DenseBins<ParticleTileDataType>;
     using index_type = ParticleBins::index_type;
     using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
 
@@ -154,12 +157,13 @@ public:
             // other species and we need to decrease their weight accordingly.
             // c1 corresponds to the minimum number of times a particle of species 1 will be paired
             // with a particle of species 2. Same for c2.
-            const index_type c1 = amrex::max(NI2/NI1,1u);
-            const index_type c2 = amrex::max(NI1/NI2,1u);
+            // index_type(1): https://github.com/AMReX-Codes/amrex/pull/3684
+            const index_type c1 = amrex::max(NI2/NI1, index_type(1));
+            const index_type c2 = amrex::max(NI1/NI2, index_type(1));
 
             // multiplier ratio to take into account unsampled pairs
             const auto multiplier_ratio = static_cast<int>(
-                (m_isSameSpecies)?(2u*max_N - 1):(max_N));
+                m_isSameSpecies ? 2*max_N - 1 : max_N);
 
 #if (defined WARPX_DIM_RZ)
             amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta];
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H
index 7b29267ec32..0b51d6b4b61 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H
@@ -22,10 +22,12 @@
 
 namespace {
     // Define shortcuts for frequently-used type names
-    using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
-    using ParticleType = WarpXParticleContainer::ParticleType;
-    using ParticleBins = amrex::DenseBins<ParticleType>;
-    using index_type = ParticleBins::index_type;
+    using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
+    using ParticleType = typename WarpXParticleContainer::ParticleType;
+    using ParticleTileType = typename WarpXParticleContainer::ParticleTileType;
+    using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType;
+    using ParticleBins = amrex::DenseBins<ParticleTileDataType>;
+    using index_type = typename ParticleBins::index_type;
 
     /**
      * \brief This function initializes the momentum of the alpha particles produced from
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H
index be3f5b2d957..52e9db8aa94 100644
--- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H
+++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/TwoProductFusionInitializeMomentum.H
@@ -24,7 +24,9 @@ namespace {
     // Define shortcuts for frequently-used type names
     using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
     using ParticleType = WarpXParticleContainer::ParticleType;
-    using ParticleBins = amrex::DenseBins<ParticleType>;
+    using ParticleTileType = WarpXParticleContainer::ParticleTileType;
+    using ParticleTileDataType = ParticleTileType::ParticleTileDataType;
+    using ParticleBins = amrex::DenseBins<ParticleTileDataType>;
     using index_type = ParticleBins::index_type;
 
     /**
diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
index dc830b477df..3928c74223d 100644
--- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
+++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H
@@ -30,13 +30,15 @@
  * \brief This functor creates particles produced from a binary collision and sets their initial
  * properties (position, momentum, weight).
  */
-class ParticleCreationFunc{
+class ParticleCreationFunc
+{
     // Define shortcuts for frequently-used type names
-    using ParticleType = WarpXParticleContainer::ParticleType;
-    using ParticleTileType = WarpXParticleContainer::ParticleTileType;
-    using ParticleBins = amrex::DenseBins<ParticleType>;
-    using index_type = ParticleBins::index_type;
-    using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
+    using ParticleType = typename WarpXParticleContainer::ParticleType;
+    using ParticleTileType = typename WarpXParticleContainer::ParticleTileType;
+    using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType;
+    using ParticleBins = amrex::DenseBins<ParticleTileDataType>;
+    using index_type = typename ParticleBins::index_type;
+    using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
 
 public:
     /**
@@ -69,12 +71,6 @@ public:
      * @param[in, out] soa_1 struct of array data of the first colliding particle species
      * @param[in, out] soa_2 struct of array data of the second colliding particle species
      * @param[out] tile_products array containing tile data of the product particles.
-     * @param[out] particle_ptr_1 pointer to data of the first colliding particle species. Is
-     *             needed to set the id of a particle to -1 in order to delete it when its weight
-     *             reaches 0.
-     * @param[out] particle_ptr_2 pointer to data of the second colliding particle species. Is
-     *             needed to set the id of a particle to -1 in order to delete it when its weight
-     *             reaches 0.
      * @param[in] m1 mass of the first colliding particle species
      * @param[in] m2 mass of the second colliding particle species
      * @param[in] products_mass array storing the mass of product particles
@@ -102,7 +98,6 @@ public:
                     const SoaData_type& soa_1, const SoaData_type& soa_2,
                     const amrex::Vector<WarpXParticleContainer*>& pc_products,
                     ParticleTileType** AMREX_RESTRICT tile_products,
-                    ParticleType* particle_ptr_1, ParticleType*  particle_ptr_2,
                     const amrex::ParticleReal& m1, const amrex::ParticleReal& m2,
                     const amrex::Vector<amrex::ParticleReal>& products_mass,
                     const index_type* AMREX_RESTRICT p_mask,
@@ -137,6 +132,8 @@ public:
 
         amrex::ParticleReal* AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
         amrex::ParticleReal* AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
+        uint64_t* AMREX_RESTRICT idcpu1 = soa_1.m_idcpu;
+        uint64_t* AMREX_RESTRICT idcpu2 = soa_2.m_idcpu;
 
         // Create necessary GPU vectors, that will be used in the kernel below
         amrex::Vector<SoaData_type> soa_products;
@@ -205,16 +202,32 @@ public:
                 amrex::Gpu::Atomic::AddNoRet(&w2[p_pair_indices_2[i]],
                                                 -p_pair_reaction_weight[i]);
 
+                // Note: Particle::atomicSetID should also be provided as a standalone helper function in AMReX
+                //       to replace the following lambda.
+                auto const atomicSetIdMinus = [] AMREX_GPU_DEVICE (uint64_t & idcpu)
+                {
+                    constexpr amrex::Long minus_one_long = -1;
+                    uint64_t tmp = 0;
+                    amrex::ParticleIDWrapper wrapper(tmp);
+                    wrapper = minus_one_long;
+#if defined(AMREX_USE_OMP)
+#pragma omp atomic write
+                    idcpu = wrapper.m_idata;
+#else
+                    auto *old_ptr = reinterpret_cast<unsigned long long*>(&idcpu);
+                    amrex::Gpu::Atomic::Exch(old_ptr, (unsigned long long) wrapper.m_idata);
+#endif
+                };
+
                 // If the colliding particle weight decreases to zero, remove particle by
                 // setting its id to -1
-                constexpr amrex::Long minus_one_long = -1;
                 if (w1[p_pair_indices_1[i]] <= amrex::ParticleReal(0.))
                 {
-                    particle_ptr_1[p_pair_indices_1[i]].atomicSetID(minus_one_long);
+                    atomicSetIdMinus(idcpu1[p_pair_indices_1[i]]);
                 }
                 if (w2[p_pair_indices_2[i]] <= amrex::ParticleReal(0.))
                 {
-                    particle_ptr_2[p_pair_indices_2[i]].atomicSetID(minus_one_long);
+                    atomicSetIdMinus(idcpu2[p_pair_indices_2[i]]);
                 }
 
                 // Initialize the product particles' momentum, using a function depending on the
@@ -294,12 +307,14 @@ private:
  * \brief This class does nothing and is used as second template parameter for binary collisions
  * that do not create particles.
  */
-class NoParticleCreationFunc{
-    using ParticleType = WarpXParticleContainer::ParticleType;
-    using ParticleTileType = WarpXParticleContainer::ParticleTileType;
-    using ParticleBins = amrex::DenseBins<ParticleType>;
-    using index_type = ParticleBins::index_type;
-    using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
+class NoParticleCreationFunc
+{
+    using ParticleType = typename WarpXParticleContainer::ParticleType;
+    using ParticleTileType = typename WarpXParticleContainer::ParticleTileType;
+    using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType;
+    using ParticleBins = amrex::DenseBins<ParticleTileDataType>;
+    using index_type = typename ParticleBins::index_type;
+    using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
 
 public:
     NoParticleCreationFunc () = default;
@@ -313,7 +328,6 @@ public:
                     const SoaData_type& /*soa_1*/, const SoaData_type& /*soa_2*/,
                     amrex::Vector<WarpXParticleContainer*>& /*pc_products*/,
                     ParticleTileType** /*tile_products*/,
-                    ParticleType* /*particle_ptr_1*/, ParticleType* /*particle_ptr_2*/,
                     const amrex::ParticleReal& /*m1*/, const amrex::ParticleReal& /*m2*/,
                     const amrex::Vector<amrex::ParticleReal>& /*products_mass*/,
                     const index_type* /*p_mask*/, const amrex::Vector<index_type>& /*products_np*/,
diff --git a/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H b/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H
index 42259512b0d..3b8f72f4b84 100644
--- a/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H
+++ b/Source/Particles/Collision/BinaryCollision/ShuffleFisherYates.H
@@ -12,7 +12,7 @@
 /* \brief Shuffle array according to Fisher-Yates algorithm.
  *        Only shuffle the part between is <= i < ie, n = ie-is.
  *        T_index shall be
- *        amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type
+ *        amrex::DenseBins<WarpXParticleContainer::ParticleTileType::ParticleTileDataType>::index_type
 */
 
 template <typename T_index>
diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H
index d0db678dfda..d0822789015 100644
--- a/Source/Particles/Deposition/ChargeDeposition.H
+++ b/Source/Particles/Deposition/ChargeDeposition.H
@@ -252,7 +252,7 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition<PIdx>& GetPositio
                                      const int n_rz_azimuthal_modes,
                                      amrex::Real* cost,
                                      const long load_balance_costs_update_algo,
-                                     const amrex::DenseBins<WarpXParticleContainer::ParticleType>& a_bins,
+                                     const amrex::DenseBins<WarpXParticleContainer::ParticleTileType::ParticleTileDataType>& a_bins,
                                      const amrex::Box& box,
                                      const amrex::Geometry& geom,
                                      const amrex::IntVect& a_tbox_max_size,
diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H
index 18df09c3b43..2252a63fd07 100644
--- a/Source/Particles/Deposition/CurrentDeposition.H
+++ b/Source/Particles/Deposition/CurrentDeposition.H
@@ -592,7 +592,7 @@ void doDepositionSharedShapeN (const GetParticlePosition<PIdx>& GetPosition,
                                int n_rz_azimuthal_modes,
                                amrex::Real* cost,
                                long load_balance_costs_update_algo,
-                               const amrex::DenseBins<WarpXParticleContainer::ParticleType>& a_bins,
+                               const amrex::DenseBins<WarpXParticleContainer::ParticleTileType::ParticleTileDataType>& a_bins,
                                const amrex::Box& box,
                                const amrex::Geometry& geom,
                                const amrex::IntVect& a_tbox_max_size)
diff --git a/Source/Particles/ElementaryProcess/QEDPairGeneration.H b/Source/Particles/ElementaryProcess/QEDPairGeneration.H
index 5abc9282d4f..ca00b56323a 100644
--- a/Source/Particles/ElementaryProcess/QEDPairGeneration.H
+++ b/Source/Particles/ElementaryProcess/QEDPairGeneration.H
@@ -167,7 +167,7 @@ public:
             p_ux, p_uy, p_uz,
             engine);
 
-        src.m_aos[i_src].id() = -1; //destroy photon after pair generation
+        amrex::ParticleIDWrapper{src.m_idcpu[i_src]} = -1;  // destroy photon after pair generation
     }
 
 private:
diff --git a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H
index 8ba5c63ad57..f3099bf997f 100644
--- a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H
+++ b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H
@@ -22,6 +22,7 @@
 #include <AMReX_GpuLaunch.H>
 #include <AMReX_GpuQualifiers.H>
 #include <AMReX_IndexType.H>
+#include <AMReX_ParticleTile.H>
 #include <AMReX_REAL.H>
 
 #include <AMReX_BaseFwd.H>
@@ -237,12 +238,11 @@ void cleanLowEnergyPhotons(
     const int old_size, const int num_added,
     const amrex::ParticleReal energy_threshold)
 {
-    auto pp = ptile.GetArrayOfStructs()().data() + old_size;
-
-    const auto& soa = ptile.GetStructOfArrays();
+    auto& soa = ptile.GetStructOfArrays();
+    auto p_idcpu = soa.GetIdCPUData().data() + old_size;
     const auto p_ux = soa.GetRealData(PIdx::ux).data() + old_size;
-    const auto  p_uy = soa.GetRealData(PIdx::uy).data() + old_size;
-    const auto  p_uz = soa.GetRealData(PIdx::uz).data() + old_size;
+    const auto p_uy = soa.GetRealData(PIdx::uy).data() + old_size;
+    const auto p_uz = soa.GetRealData(PIdx::uz).data() + old_size;
 
     //The square of the energy threshold
     const auto energy_threshold2 = std::max(
@@ -251,8 +251,6 @@ void cleanLowEnergyPhotons(
 
     amrex::ParallelFor(num_added, [=] AMREX_GPU_DEVICE (int ip) noexcept
     {
-        auto& p = pp[ip];
-
         const auto ux = p_ux[ip];
         const auto uy = p_uy[ip];
         const auto uz = p_uz[ip];
@@ -262,8 +260,8 @@ void cleanLowEnergyPhotons(
         constexpr amrex::ParticleReal me_c = PhysConst::m_e*PhysConst::c;
         const auto phot_energy2 = (ux*ux + uy*uy + uz*uz)*me_c*me_c;
 
-        if (phot_energy2 < energy_threshold2){
-            p.id() = - 1;
+        if (phot_energy2 < energy_threshold2) {
+            amrex::ParticleIDWrapper{p_idcpu[ip]} = -1;
         }
     });
 }
diff --git a/Source/Particles/LaserParticleContainer.H b/Source/Particles/LaserParticleContainer.H
index e6fa308431c..fac94ff20a3 100644
--- a/Source/Particles/LaserParticleContainer.H
+++ b/Source/Particles/LaserParticleContainer.H
@@ -56,10 +56,9 @@ public:
      * \brief Method to initialize runtime attributes. Does nothing for LaserParticleContainer.
      */
     void DefaultInitializeRuntimeAttributes (
-                        amrex::ParticleTile<amrex::Particle<NStructReal, NStructInt>,
-                                            NArrayReal, NArrayInt, amrex::PinnedArenaAllocator>& /*pinned_tile*/,
-                        const int /*n_external_attr_real*/,
-                        const int /*n_external_attr_int*/) final {}
+        typename ContainerLike<amrex::PinnedArenaAllocator>::ParticleTileType& /*pinned_tile*/,
+        int /*n_external_attr_real*/,
+        int /*n_external_attr_int*/) final {}
 
     void ReadHeader (std::istream& is) final;
 
diff --git a/Source/Particles/NamedComponentParticleContainer.H b/Source/Particles/NamedComponentParticleContainer.H
index 3be0886425d..e7a7a20fad5 100644
--- a/Source/Particles/NamedComponentParticleContainer.H
+++ b/Source/Particles/NamedComponentParticleContainer.H
@@ -18,24 +18,39 @@
 #include <utility>
 
 
-/** Particle Attributes stored in amrex::ParticleContainer's struct of array
+/** Real Particle Attributes stored in amrex::ParticleContainer's struct of array
  */
 struct PIdx
 {
     enum {
-        w = 0,      ///< weight
+#if !defined (WARPX_DIM_1D_Z)
+        x,
+#endif
+#if defined (WARPX_DIM_3D)
+        y,
+#endif
+        z,
+        w,      ///< weight
         ux, uy, uz,
 #ifdef WARPX_DIM_RZ
         theta,      ///< RZ needs all three position components
 #endif
-        nattribs    ///< number of attributes
+        nattribs    ///< number of compile-time attributes
+    };
+};
+
+/** Integer Particle Attributes stored in amrex::ParticleContainer's struct of array
+ */
+struct PIdxInt
+{
+    enum {
+        nattribs    ///< number of compile-time attributes
     };
 };
 
 /** Particle Container class that allows to add/access particle components
  *  with a name (string) instead of doing so with an integer index.
- *  (The "components" are all the particle quantities - except those
- *  that are stored in an AoS by amrex, i.e. the particle positions and ID)
+ *  (The "components" are all the particle amrex::Real quantities.)
  *
  *  This is done by storing maps that give the index of the component
  *  that corresponds to a given string.
@@ -45,11 +60,11 @@ struct PIdx
  */
 template <template<class> class T_Allocator=amrex::DefaultAllocator>
 class NamedComponentParticleContainer :
-public amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>
+public amrex::ParticleContainerPureSoA<PIdx::nattribs, 0, T_Allocator>
 {
 public:
     /** Construct an empty NamedComponentParticleContainer **/
-    NamedComponentParticleContainer () : amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>() {}
+    NamedComponentParticleContainer () : amrex::ParticleContainerPureSoA<PIdx::nattribs, 0, T_Allocator>() {}
 
     /** Construct a NamedComponentParticleContainer from an AmrParGDB object
      *
@@ -61,8 +76,15 @@ public:
      * AMR hierarchy. Usually, this is generated by an AmrCore or AmrLevel object.
      */
     NamedComponentParticleContainer (amrex::AmrParGDB* amr_pgdb)
-    : amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>(amr_pgdb) {
+    : amrex::ParticleContainerPureSoA<PIdx::nattribs, 0, T_Allocator>(amr_pgdb) {
         // build up the map of string names to particle component numbers
+#if !defined (WARPX_DIM_1D_Z)
+        particle_comps["x"]  = PIdx::x;
+#endif
+#if defined (WARPX_DIM_3D)
+        particle_comps["y"]  = PIdx::y;
+#endif
+        particle_comps["z"]  = PIdx::z;
         particle_comps["w"]  = PIdx::w;
         particle_comps["ux"] = PIdx::ux;
         particle_comps["uy"] = PIdx::uy;
@@ -85,12 +107,12 @@ public:
      * @param p_ricomps name-to-index map for run-time integer components
      */
     NamedComponentParticleContainer(
-        amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator> && pc,
+        amrex::ParticleContainerPureSoA<PIdx::nattribs, 0, T_Allocator> && pc,
         std::map<std::string, int> p_comps,
         std::map<std::string, int> p_icomps,
         std::map<std::string, int> p_rcomps,
         std::map<std::string, int> p_ricomps)
-    : amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>(std::move(pc)),
+    : amrex::ParticleContainerPureSoA<PIdx::nattribs, 0, T_Allocator>(std::move(pc)),
     particle_comps(std::move(p_comps)),
     particle_icomps(std::move(p_icomps)),
     particle_runtime_comps(std::move(p_rcomps)),
@@ -118,7 +140,7 @@ public:
     NamedComponentParticleContainer<NewAllocator>
     make_alike () const {
         auto tmp = NamedComponentParticleContainer<NewAllocator>(
-            amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::template make_alike<NewAllocator>(),
+            amrex::ParticleContainerPureSoA<PIdx::nattribs, 0, T_Allocator>::template make_alike<NewAllocator>(),
             particle_comps,
             particle_icomps,
             particle_runtime_comps,
@@ -127,10 +149,10 @@ public:
         return tmp;
     }
 
-    using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::NumRealComps;
-    using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::NumIntComps;
-    using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::AddRealComp;
-    using amrex::ParticleContainer<0,0,PIdx::nattribs,0,T_Allocator>::AddIntComp;
+    using amrex::ParticleContainerPureSoA<PIdx::nattribs,0,T_Allocator>::NumRealComps;
+    using amrex::ParticleContainerPureSoA<PIdx::nattribs,0,T_Allocator>::NumIntComps;
+    using amrex::ParticleContainerPureSoA<PIdx::nattribs,0,T_Allocator>::AddRealComp;
+    using amrex::ParticleContainerPureSoA<PIdx::nattribs,0,T_Allocator>::AddIntComp;
 
     /** Allocate a new run-time real component
      *
diff --git a/Source/Particles/ParticleBoundaryBuffer.cpp b/Source/Particles/ParticleBoundaryBuffer.cpp
index 54c4396379d..88304bd8a9c 100644
--- a/Source/Particles/ParticleBoundaryBuffer.cpp
+++ b/Source/Particles/ParticleBoundaryBuffer.cpp
@@ -50,7 +50,7 @@ struct CopyAndTimestamp {
     void operator() (const DstData& dst, const SrcData& src,
                      int src_i, int dst_i) const noexcept
     {
-        dst.m_aos[dst_i] = src.m_aos[src_i];
+        dst.m_idcpu[dst_i] = src.m_idcpu[src_i];
         for (int j = 0; j < SrcData::NAR; ++j) {
             dst.m_rdata[j][dst_i] = src.m_rdata[j][src_i];
         }
@@ -222,7 +222,7 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc,
 {
     WARPX_PROFILE("ParticleBoundaryBuffer::gatherParticles");
 
-    using PIter = amrex::ParConstIter<0,0,PIdx::nattribs>;
+    using PIter = amrex::ParConstIterSoA<PIdx::nattribs, 0>;
     const auto& warpx_instance = WarpX::GetInstance();
     const amrex::Geometry& geom = warpx_instance.Geom(0);
     auto plo = geom.ProbLoArray();
diff --git a/Source/Particles/ParticleCreation/SmartCopy.H b/Source/Particles/ParticleCreation/SmartCopy.H
index 2c04baa18bb..6a6ceb3d290 100644
--- a/Source/Particles/ParticleCreation/SmartCopy.H
+++ b/Source/Particles/ParticleCreation/SmartCopy.H
@@ -26,7 +26,7 @@
  * type. Second, if a given component name is found in both the src
  * and the dst, then the src value is copied.
  *
- * Particle structs - positions and id numbers  - are always copied.
+ * Particle positions and id numbers are always copied.
  *
  * You don't create this directly - use the SmartCopyFactory object below.
  */
@@ -48,9 +48,6 @@ struct SmartCopy
     void operator() (DstData& dst, const SrcData& src, int i_src, int i_dst,
                      amrex::RandomEngine const& engine) const noexcept
     {
-        // the particle struct is always copied over
-        dst.m_aos[i_dst] = src.m_aos[i_src];
-
         // initialize the real components
         for (int j = 0; j < DstData::NAR; ++j) {
             dst.m_rdata[j][i_dst] = initializeRealValue(m_policy_real[j], engine);
diff --git a/Source/Particles/ParticleCreation/SmartCreate.H b/Source/Particles/ParticleCreation/SmartCreate.H
index 67d7767a5d3..7fb854ee083 100644
--- a/Source/Particles/ParticleCreation/SmartCreate.H
+++ b/Source/Particles/ParticleCreation/SmartCreate.H
@@ -14,6 +14,8 @@
 #include <AMReX_AmrCore.H>
 #include <AMReX_GpuContainers.H>
 #include <AMReX_ParallelDescriptor.H>
+#include <AMReX_Particle.H>
+#include <AMReX_ParticleTile.H>
 
 /**
  * \brief This is a functor for performing a "smart create" that works
@@ -47,23 +49,23 @@ struct SmartCreate
         const int id = 0) const noexcept
     {
 #if defined(WARPX_DIM_3D)
-        prt.m_aos[i_prt].pos(0) = x;
-        prt.m_aos[i_prt].pos(1) = y;
-        prt.m_aos[i_prt].pos(2) = z;
+        prt.m_rdata[PIdx::x][i_prt] = x;
+        prt.m_rdata[PIdx::y][i_prt] = y;
+        prt.m_rdata[PIdx::z][i_prt] = z;
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-        prt.m_aos[i_prt].pos(0) = x;
-        prt.m_aos[i_prt].pos(1) = z;
+        prt.m_rdata[PIdx::x][i_prt] = x;
+        prt.m_rdata[PIdx::z][i_prt] = z;
         amrex::ignore_unused(y);
 #else
-        prt.m_aos[i_prt].pos(0) = z;
+        prt.m_rdata[PIdx::z][i_prt] = z;
         amrex::ignore_unused(x,y);
 #endif
 
-        prt.m_aos[i_prt].cpu() = cpu;
-        prt.m_aos[i_prt].id() = id;
+        amrex::ParticleIDWrapper{prt.m_idcpu[i_prt]} = id;
+        amrex::ParticleCPUWrapper{prt.m_idcpu[i_prt]} = cpu;
 
-        // initialize the real components
-        for (int j = 0; j < PartData::NAR; ++j) {
+        // initialize the real components after position
+        for (int j = AMREX_SPACEDIM; j < PartData::NAR; ++j) {
             prt.m_rdata[j][i_prt] = initializeRealValue(m_policy_real[j], engine);
         }
         for (int j = 0; j < prt.m_num_runtime_real; ++j) {
diff --git a/Source/Particles/ParticleCreation/SmartUtils.H b/Source/Particles/ParticleCreation/SmartUtils.H
index 732a12bb729..6b3d396900d 100644
--- a/Source/Particles/ParticleCreation/SmartUtils.H
+++ b/Source/Particles/ParticleCreation/SmartUtils.H
@@ -17,6 +17,7 @@
 #include <AMReX_GpuQualifiers.H>
 #include <AMReX_INT.H>
 #include <AMReX_ParallelDescriptor.H>
+#include <AMReX_Particle.H>
 
 #include <map>
 #include <string>
@@ -60,12 +61,12 @@ void setNewParticleIDs (PTile& ptile, int old_size, int num_added)
     }
 
     const int cpuid = amrex::ParallelDescriptor::MyProc();
-    auto pp = ptile.GetArrayOfStructs()().data() + old_size;
+    auto ptd = ptile.getParticleTileData();
     amrex::ParallelFor(num_added, [=] AMREX_GPU_DEVICE (int ip) noexcept
     {
-        auto& p = pp[ip];
-        p.id() = pid+ip;
-        p.cpu() = cpuid;
+        auto const new_id = ip + old_size;
+        amrex::ParticleIDWrapper{ptd.m_idcpu[new_id]} = pid+ip;
+        amrex::ParticleCPUWrapper{ptd.m_idcpu[new_id]} = cpuid;
     });
 }
 
diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H
index a12ae75f629..edf91a84526 100644
--- a/Source/Particles/PhysicalParticleContainer.H
+++ b/Source/Particles/PhysicalParticleContainer.H
@@ -268,11 +268,9 @@ public:
      * @param[in] engine the random engine, used in initialization of QED optical depths
      */
     void DefaultInitializeRuntimeAttributes (
-                    amrex::ParticleTile<amrex::Particle<NStructReal, NStructInt>,
-                                        NArrayReal, NArrayInt,
-                                        amrex::PinnedArenaAllocator>& pinned_tile,
-                    int n_external_attr_real,
-                    int n_external_attr_int) final;
+        typename ContainerLike<amrex::PinnedArenaAllocator>::ParticleTileType& pinned_tile,
+        int n_external_attr_real,
+        int n_external_attr_int) final;
 
 /**
  * \brief Apply NCI Godfrey filter to all components of E and B before gather
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 929c3c26649..bef78d714bd 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -198,8 +198,8 @@ namespace
      * and avoid any possible undefined behavior before the next call to redistribute) and sets
      * the particle id to -1 so that it can be effectively deleted.
      *
-     * \param p particle aos data
-     * \param pa particle soa data
+     * \param idcpu particle id soa data
+     * \param pa particle real soa data
      * \param ip index for soa data
      * \param do_field_ionization whether species has ionization
      * \param pi ionization level data
@@ -210,20 +210,21 @@ namespace
      */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     void ZeroInitializeAndSetNegativeID (
-        ParticleType& p, const GpuArray<ParticleReal*,PIdx::nattribs>& pa, long& ip,
+        uint64_t * AMREX_RESTRICT idcpu,
+        const GpuArray<ParticleReal*,PIdx::nattribs>& pa, long& ip,
         const bool& do_field_ionization, int* pi
 #ifdef WARPX_QED
-        ,const bool& has_quantum_sync, amrex::ParticleReal* p_optical_depth_QSR
-        ,const bool& has_breit_wheeler, amrex::ParticleReal* p_optical_depth_BW
+        ,const bool& has_quantum_sync, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_QSR
+        ,const bool& has_breit_wheeler, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_BW
 #endif
         ) noexcept
     {
-        p.pos(0) = 0._rt;
+        pa[PIdx::z][ip] = 0._rt;
 #if (AMREX_SPACEDIM >= 2)
-        p.pos(1) = 0._rt;
+        pa[PIdx::x][ip] = 0._rt;
 #endif
 #if defined(WARPX_DIM_3D)
-        p.pos(2) = 0._rt;
+        pa[PIdx::y][ip] = 0._rt;
 #endif
         pa[PIdx::w ][ip] = 0._rt;
         pa[PIdx::ux][ip] = 0._rt;
@@ -238,7 +239,7 @@ namespace
         if (has_breit_wheeler) {p_optical_depth_BW[ip] = 0._rt;}
 #endif
 
-        p.id() = -1;
+        amrex::ParticleIDWrapper{idcpu[ip]} = -1;
     }
 }
 
@@ -780,11 +781,9 @@ PhysicalParticleContainer::AddPlasmaFromFile(PlasmaInjector & plasma_injector,
 
 void
 PhysicalParticleContainer::DefaultInitializeRuntimeAttributes (
-                    amrex::ParticleTile<amrex::Particle<NStructReal, NStructInt>,
-                                        NArrayReal, NArrayInt,
-                                        amrex::PinnedArenaAllocator>& pinned_tile,
-                    const int n_external_attr_real,
-                    const int n_external_attr_int)
+    typename ContainerLike<amrex::PinnedArenaAllocator>::ParticleTileType& pinned_tile,
+    int n_external_attr_real,
+    int n_external_attr_int)
 {
     ParticleCreation::DefaultInitializeRuntimeAttributes(pinned_tile,
                                        n_external_attr_real, n_external_attr_int,
@@ -1084,7 +1083,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
         const int max_new_particles = Scan::ExclusiveSum(counts.size(), counts.data(), offset.data());
 
         // Update NextID to include particles created in this function
-        Long pid;
+        int pid;
 #ifdef AMREX_USE_OMP
 #pragma omp critical (add_plasma_nextid)
 #endif
@@ -1093,7 +1092,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
             ParticleType::NextID(pid+max_new_particles);
         }
         WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            static_cast<Long>(pid + max_new_particles) < LastParticleID,
+            static_cast<amrex::Long>(pid) + static_cast<amrex::Long>(max_new_particles) < LongParticleIds::LastParticleID,
             "ERROR: overflow on particle id numbers");
 
         const int cpuid = ParallelDescriptor::MyProc();
@@ -1104,16 +1103,16 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
             DefineAndReturnParticleTile(lev, grid_id, tile_id);
         }
 
-        auto old_size = particle_tile.GetArrayOfStructs().size();
+        auto old_size = particle_tile.size();
         auto new_size = old_size + max_new_particles;
         particle_tile.resize(new_size);
 
-        ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size;
         auto& soa = particle_tile.GetStructOfArrays();
         GpuArray<ParticleReal*,PIdx::nattribs> pa;
         for (int ia = 0; ia < PIdx::nattribs; ++ia) {
             pa[ia] = soa.GetRealData(ia).data() + old_size;
         }
+        uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size;
         // user-defined integer and real attributes
         const auto n_user_int_attribs = static_cast<int>(m_user_int_attribs.size());
         const auto n_user_real_attribs = static_cast<int>(m_user_real_attribs.size());
@@ -1226,9 +1225,8 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
             for (int i_part = 0; i_part < pcounts[index]; ++i_part)
             {
                 long ip = poffset[index] + i_part;
-                ParticleType& p = pp[ip];
-                p.id() = pid+ip;
-                p.cpu() = cpuid;
+                amrex::ParticleIDWrapper{pa_idcpu[ip]} = pid+ip;
+                amrex::ParticleCPUWrapper{pa_idcpu[ip]} = cpuid;
                 const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ?
                   // In the refined injection region: use refinement ratio `lrrfac`
                   inj_pos->getPositionUnitBox(i_part, lrrfac, engine) :
@@ -1238,7 +1236,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
 
 #if defined(WARPX_DIM_3D)
                 if (!tile_realbox.contains(XDim3{pos.x,pos.y,pos.z})) {
-                    ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+                    ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi
 #ifdef WARPX_QED
                                                    ,loc_has_quantum_sync, p_optical_depth_QSR
                                                    ,loc_has_breit_wheeler, p_optical_depth_BW
@@ -1249,7 +1247,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
                 amrex::ignore_unused(k);
                 if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) {
-                    ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+                    ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi
 #ifdef WARPX_QED
                                                    ,loc_has_quantum_sync, p_optical_depth_QSR
                                                    ,loc_has_breit_wheeler, p_optical_depth_BW
@@ -1260,7 +1258,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
 #else
                 amrex::ignore_unused(j,k);
                 if (!tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt})) {
-                    ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+                    ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi
 #ifdef WARPX_QED
                                                    ,loc_has_quantum_sync, p_optical_depth_QSR
                                                    ,loc_has_breit_wheeler, p_optical_depth_BW
@@ -1299,7 +1297,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
                     const Real z0 = applyBallisticCorrection(pos, inj_mom, gamma_boost,
                                                              beta_boost, t);
                     if (!inj_pos->insideBounds(xb, yb, z0)) {
-                        ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+                        ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi
 #ifdef WARPX_QED
                                                    ,loc_has_quantum_sync, p_optical_depth_QSR
                                                    ,loc_has_breit_wheeler, p_optical_depth_BW
@@ -1313,7 +1311,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
 
                     // Remove particle if density below threshold
                     if ( dens < density_min ){
-                        ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+                        ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi
 #ifdef WARPX_QED
                                                    ,loc_has_quantum_sync, p_optical_depth_QSR
                                                    ,loc_has_breit_wheeler, p_optical_depth_BW
@@ -1331,7 +1329,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
                     // If the particle is not within the lab-frame zmin, zmax, etc.
                     // go to the next generated particle.
                     if (!inj_pos->insideBounds(xb, yb, z0_lab)) {
-                        ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+                        ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi
 #ifdef WARPX_QED
                                                    ,loc_has_quantum_sync, p_optical_depth_QSR
                                                    ,loc_has_breit_wheeler, p_optical_depth_BW
@@ -1343,7 +1341,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
                     dens = inj_rho->getDensity(pos.x, pos.y, z0_lab);
                     // Remove particle if density below threshold
                     if ( dens < density_min ){
-                        ZeroInitializeAndSetNegativeID(p, pa, ip, loc_do_field_ionization, pi
+                        ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi
 #ifdef WARPX_QED
                                                    ,loc_has_quantum_sync, p_optical_depth_QSR
                                                    ,loc_has_breit_wheeler, p_optical_depth_BW
@@ -1410,17 +1408,17 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int
                 pa[PIdx::uz][ip] = u.z;
 
 #if defined(WARPX_DIM_3D)
-                p.pos(0) = pos.x;
-                p.pos(1) = pos.y;
-                p.pos(2) = pos.z;
+                pa[PIdx::x][ip] = pos.x;
+                pa[PIdx::y][ip] = pos.y;
+                pa[PIdx::z][ip] = pos.z;
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
 #ifdef WARPX_DIM_RZ
                 pa[PIdx::theta][ip] = theta;
 #endif
-                p.pos(0) = xb;
-                p.pos(1) = pos.z;
+                pa[PIdx::x][ip] = xb;
+                pa[PIdx::z][ip] = pos.z;
 #else
-                p.pos(0) = pos.z;
+                pa[PIdx::z][ip] = pos.z;
 #endif
             }
         });
@@ -1645,7 +1643,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
         const int max_new_particles = Scan::ExclusiveSum(counts.size(), counts.data(), offset.data());
 
         // Update NextID to include particles created in this function
-        Long pid;
+        int pid;
 #ifdef AMREX_USE_OMP
 #pragma omp critical (add_plasma_nextid)
 #endif
@@ -1654,23 +1652,23 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
             ParticleType::NextID(pid+max_new_particles);
         }
         WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            static_cast<Long>(pid + max_new_particles) < LastParticleID,
+            static_cast<amrex::Long>(pid) + static_cast<amrex::Long>(max_new_particles) < LongParticleIds::LastParticleID,
             "overflow on particle id numbers");
 
         const int cpuid = ParallelDescriptor::MyProc();
 
         auto& particle_tile = tmp_pc.DefineAndReturnParticleTile(0, grid_id, tile_id);
 
-        auto old_size = particle_tile.GetArrayOfStructs().size();
+        auto old_size = particle_tile.size();
         auto new_size = old_size + max_new_particles;
         particle_tile.resize(new_size);
 
-        ParticleType* pp = particle_tile.GetArrayOfStructs()().data() + old_size;
         auto& soa = particle_tile.GetStructOfArrays();
         GpuArray<ParticleReal*,PIdx::nattribs> pa;
         for (int ia = 0; ia < PIdx::nattribs; ++ia) {
             pa[ia] = soa.GetRealData(ia).data() + old_size;
         }
+        uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size;
 
         // user-defined integer and real attributes
         const auto n_user_int_attribs = static_cast<int>(m_user_int_attribs.size());
@@ -1768,9 +1766,8 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
             for (int i_part = 0; i_part < pcounts[index]; ++i_part)
             {
                 const long ip = poffset[index] + i_part;
-                ParticleType& p = pp[ip];
-                p.id() = pid+ip;
-                p.cpu() = cpuid;
+                amrex::ParticleIDWrapper{pa_idcpu[ip]} = pid+ip;
+                amrex::ParticleCPUWrapper{pa_idcpu[ip]} = cpuid;
 
                 // This assumes the flux_pos is of type InjectorPositionRandomPlane
                 const XDim3 r = (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) ?
@@ -1795,19 +1792,19 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
                 // the particles will be within the domain.
 #if defined(WARPX_DIM_3D)
                 if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.x,ppos.y,ppos.z})) {
-                    p.id() = -1;
+                    amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1;
                     continue;
                 }
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
                 amrex::ignore_unused(k);
                 if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.x,ppos.z,0.0_prt})) {
-                    p.id() = -1;
+                    amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1;
                     continue;
                 }
 #else
                 amrex::ignore_unused(j,k);
                 if (!ParticleUtils::containsInclusive(tile_realbox, XDim3{ppos.z,0.0_prt,0.0_prt})) {
-                    p.id() = -1;
+                    amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1;
                     continue;
                 }
 #endif
@@ -1815,7 +1812,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
                 // If the particle's initial position is not within or on the species's
                 // xmin, xmax, ymin, ymax, zmin, zmax, go to the next generated particle.
                 if (!flux_pos->insideBoundsInclusive(ppos.x, ppos.y, ppos.z)) {
-                    p.id() = -1;
+                    amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1;
                     continue;
                 }
 
@@ -1849,7 +1846,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
                 Real flux = inj_flux->getFlux(ppos.x, ppos.y, ppos.z, t);
                 // Remove particle if flux is negative or 0
                 if ( flux <=0 ){
-                    p.id() = -1;
+                    amrex::ParticleIDWrapper{pa_idcpu[ip]} = -1;
                     continue;
                 }
 
@@ -1908,18 +1905,18 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector,
                 UpdatePosition(ppos.x, ppos.y, ppos.z, pu.x, pu.y, pu.z, t_fract);
 
 #if defined(WARPX_DIM_3D)
-                p.pos(0) = ppos.x;
-                p.pos(1) = ppos.y;
-                p.pos(2) = ppos.z;
+                pa[PIdx::x][ip] = ppos.x;
+                pa[PIdx::y][ip] = ppos.y;
+                pa[PIdx::z][ip] = ppos.z;
 #elif defined(WARPX_DIM_RZ)
                 pa[PIdx::theta][ip] = std::atan2(ppos.y, ppos.x);
-                p.pos(0) = std::sqrt(ppos.x*ppos.x + ppos.y*ppos.y);
-                p.pos(1) = ppos.z;
+                pa[PIdx::x][ip] = std::sqrt(ppos.x*ppos.x + ppos.y*ppos.y);
+                pa[PIdx::z][ip] = ppos.z;
 #elif defined(WARPX_DIM_XZ)
-                p.pos(0) = ppos.x;
-                p.pos(1) = ppos.z;
+                pa[PIdx::x][ip] = ppos.x;
+                pa[PIdx::z][ip] = ppos.z;
 #else
-                p.pos(0) = ppos.z;
+                pa[PIdx::z][ip] = ppos.z;
 #endif
             }
         });
@@ -2342,20 +2339,22 @@ PhysicalParticleContainer::SplitParticles (int lev)
             split_offset[1] /= ppc_nd[1];
             split_offset[2] /= ppc_nd[2];
         }
-        // particle Array Of Structs data
-        auto& particles = pti.GetArrayOfStructs();
         // particle Struct Of Arrays data
         auto& attribs = pti.GetAttribs();
         auto& wp  = attribs[PIdx::w ];
         auto& uxp = attribs[PIdx::ux];
         auto& uyp = attribs[PIdx::uy];
         auto& uzp = attribs[PIdx::uz];
+
+        ParticleTileType& ptile = ParticlesAt(lev, pti);
+        auto& soa = ptile.GetStructOfArrays();
+        uint64_t * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data();
+
         const long np = pti.numParticles();
         for(int i=0; i<np; i++){
             ParticleReal xp, yp, zp;
             GetPosition(i, xp, yp, zp);
-            auto& p = particles[i];
-            if (p.id() == DoSplitParticleID){
+            if (idcpu[i] == LongParticleIds::DoSplitParticleID){
                 // If particle is tagged, split it and put the
                 // split particles in local arrays psplit_x etc.
                 np_split_to_add += np_split;
@@ -2460,7 +2459,7 @@ PhysicalParticleContainer::SplitParticles (int lev)
                 }
 #endif
                 // invalidate the particle
-                p.id() = -p.id();
+                amrex::ParticleIDWrapper{idcpu[i]} = -1;
             }
         }
     }
@@ -2483,10 +2482,16 @@ PhysicalParticleContainer::SplitParticles (int lev)
     amrex::Vector<amrex::Vector<int>> attr_int;
     pctmp_split.AddNParticles(lev,
                               np_split_to_add,
-                              xp,  yp,  zp, uxp, uyp, uzp,
-                              1, attr,
+                              xp,
+                              yp,
+                              zp,
+                              uxp,
+                              uyp,
+                              uzp,
+                              1,
+                              attr,
                               0, attr_int,
-                              1, NoSplitParticleID);
+                              1, LongParticleIds::NoSplitParticleID);
     // Copy particles from tmp to current particle container
     constexpr bool local_flag = true;
     addParticles(pctmp_split,local_flag);
diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H
index e4477a2a60d..44641557756 100644
--- a/Source/Particles/Pusher/GetAndSetPosition.H
+++ b/Source/Particles/Pusher/GetAndSetPosition.H
@@ -30,24 +30,26 @@ void get_particle_position (const WarpXParticleContainer::SuperParticleType& p,
                             amrex::ParticleReal& y,
                             amrex::ParticleReal& z) noexcept
 {
-#ifdef WARPX_DIM_RZ
-    const amrex::ParticleReal theta = p.rdata(T_PIdx::theta);
-    const amrex::ParticleReal r = p.pos(0);
+    using namespace amrex::literals;
+
+#if defined(WARPX_DIM_RZ)
+    amrex::ParticleReal const theta = p.rdata(T_PIdx::theta);
+    amrex::ParticleReal const r = p.pos(T_PIdx::x);
     x = r*std::cos(theta);
     y = r*std::sin(theta);
-    z = p.pos(1);
-#elif WARPX_DIM_3D
-    x = p.pos(0);
-    y = p.pos(1);
-    z = p.pos(2);
-#elif WARPX_DIM_XZ
-    x = p.pos(0);
-    y = amrex::ParticleReal(0.0);
-    z = p.pos(1);
+    z = p.pos(PIdx::z);
+#elif defined(WARPX_DIM_3D)
+    x = p.pos(PIdx::x);
+    y = p.pos(PIdx::y);
+    z = p.pos(PIdx::z);
+#elif defined(WARPX_DIM_XZ)
+    x = p.pos(PIdx::x);
+    y = 0_prt;
+    z = p.pos(PIdx::z);
 #else
-    x = amrex::ParticleReal(0.0);
-    y = amrex::ParticleReal(0.0);
-    z = p.pos(0);
+    x = 0_prt;
+    y = 0_prt;
+    z = p.pos(PIdx::z);
 #endif
 }
 
@@ -59,10 +61,19 @@ void get_particle_position (const WarpXParticleContainer::SuperParticleType& p,
 template<typename T_PIdx = PIdx>
 struct GetParticlePosition
 {
-    using PType = WarpXParticleContainer::ParticleType;
     using RType = amrex::ParticleReal;
 
-    const PType* AMREX_RESTRICT m_structs = nullptr;
+#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
+    const RType* AMREX_RESTRICT m_x = nullptr;
+    const RType* AMREX_RESTRICT m_z = nullptr;
+#elif defined(WARPX_DIM_3D)
+    const RType* AMREX_RESTRICT m_x = nullptr;
+    const RType* AMREX_RESTRICT m_y = nullptr;
+    const RType* AMREX_RESTRICT m_z = nullptr;
+#elif defined(WARPX_DIM_1D_Z)
+    const RType* AMREX_RESTRICT m_z = nullptr;
+#endif
+
 #if defined(WARPX_DIM_RZ)
     const RType* m_theta = nullptr;
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
@@ -84,10 +95,19 @@ struct GetParticlePosition
     template <typename ptiType>
     GetParticlePosition (const ptiType& a_pti, long a_offset = 0) noexcept
     {
-        const auto& aos = a_pti.GetArrayOfStructs();
-        m_structs = aos().dataPtr() + a_offset;
-#if defined(WARPX_DIM_RZ)
         const auto& soa = a_pti.GetStructOfArrays();
+
+#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
+        m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset;
+        m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset;
+#elif defined(WARPX_DIM_3D)
+        m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset;
+        m_y = soa.GetRealData(PIdx::y).dataPtr() + a_offset;
+        m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset;
+#elif defined(WARPX_DIM_1D_Z)
+        m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset;
+#endif
+#if defined(WARPX_DIM_RZ)
         m_theta = soa.GetRealData(T_PIdx::theta).dataPtr() + a_offset;
 #endif
     }
@@ -98,24 +118,23 @@ struct GetParticlePosition
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     void operator() (const long i, RType& x, RType& y, RType& z) const noexcept
     {
-        const PType& p = m_structs[i];
 #ifdef WARPX_DIM_RZ
-        const RType r = p.pos(0);
+        RType const r = m_x[i];
         x = r*std::cos(m_theta[i]);
         y = r*std::sin(m_theta[i]);
-        z = p.pos(1);
+        z = m_z[i];
 #elif WARPX_DIM_3D
-        x = p.pos(0);
-        y = p.pos(1);
-        z = p.pos(2);
+        x = m_x[i];
+        y = m_y[i];
+        z = m_z[i];
 #elif WARPX_DIM_XZ
-        x = p.pos(0);
+        x = m_x[i];
         y = m_y_default;
-        z = p.pos(1);
+        z = m_z[i];
 #else
         x = m_x_default;
         y = m_y_default;
-        z = p.pos(0);
+        z = m_z[i];
 #endif
     }
 
@@ -127,23 +146,22 @@ struct GetParticlePosition
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     void AsStored (const long i, RType& x, RType& y, RType& z) const noexcept
     {
-        const PType& p = m_structs[i];
 #ifdef WARPX_DIM_RZ
-        x = p.pos(0);
+        x = m_x[i];
         y = m_theta[i];
-        z = p.pos(1);
+        z = m_z[i];
 #elif WARPX_DIM_3D
-        x = p.pos(0);
-        y = p.pos(1);
-        z = p.pos(2);
+        x = m_x[i];
+        y = m_y[i];
+        z = m_z[i];
 #elif WARPX_DIM_XZ
-        x = p.pos(0);
+        x = m_x[i];
         y = m_y_default;
-        z = p.pos(1);
+        z = m_z[i];
 #else
         x = m_x_default;
         y = m_y_default;
-        z = p.pos(0);
+        z = m_z[i];
 #endif
     }
 };
@@ -158,10 +176,18 @@ struct GetParticlePosition
 template<typename T_PIdx = PIdx>
 struct SetParticlePosition
 {
-    using PType = WarpXParticleContainer::ParticleType;
     using RType = amrex::ParticleReal;
 
-    PType* AMREX_RESTRICT m_structs;
+#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
+    RType* AMREX_RESTRICT m_x;
+    RType* AMREX_RESTRICT m_z;
+#elif defined(WARPX_DIM_3D)
+    RType* AMREX_RESTRICT m_x;
+    RType* AMREX_RESTRICT m_y;
+    RType* AMREX_RESTRICT m_z;
+#elif defined(WARPX_DIM_1D_Z)
+    RType* AMREX_RESTRICT m_z;
+#endif
 #if defined(WARPX_DIM_RZ)
     RType* AMREX_RESTRICT m_theta;
 #endif
@@ -169,10 +195,18 @@ struct SetParticlePosition
     template <typename ptiType>
     SetParticlePosition (const ptiType& a_pti, long a_offset = 0) noexcept
     {
-        auto& aos = a_pti.GetArrayOfStructs();
-        m_structs = aos().dataPtr() + a_offset;
-#if defined(WARPX_DIM_RZ)
         auto& soa = a_pti.GetStructOfArrays();
+#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ)
+        m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset;
+        m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset;
+#elif defined(WARPX_DIM_3D)
+        m_x = soa.GetRealData(PIdx::x).dataPtr() + a_offset;
+        m_y = soa.GetRealData(PIdx::y).dataPtr() + a_offset;
+        m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset;
+#elif defined(WARPX_DIM_1D_Z)
+        m_z = soa.GetRealData(PIdx::z).dataPtr() + a_offset;
+#endif
+#if defined(WARPX_DIM_RZ)
         m_theta = soa.GetRealData(T_PIdx::theta).dataPtr() + a_offset;
 #endif
     }
@@ -190,17 +224,17 @@ struct SetParticlePosition
 #endif
 #ifdef WARPX_DIM_RZ
         m_theta[i] = std::atan2(y, x);
-        m_structs[i].pos(0) = std::sqrt(x*x + y*y);
-        m_structs[i].pos(1) = z;
+        m_x[i] = std::sqrt(x*x + y*y);
+        m_z[i] = z;
 #elif WARPX_DIM_3D
-        m_structs[i].pos(0) = x;
-        m_structs[i].pos(1) = y;
-        m_structs[i].pos(2) = z;
+        m_x[i] = x;
+        m_y[i] = y;
+        m_z[i] = z;
 #elif WARPX_DIM_XZ
-        m_structs[i].pos(0) = x;
-        m_structs[i].pos(1) = z;
+        m_x[i] = x;
+        m_z[i] = z;
 #else
-        m_structs[i].pos(0) = z;
+        m_z[i] = z;
 #endif
     }
 
@@ -218,18 +252,18 @@ struct SetParticlePosition
         amrex::ignore_unused(x,y);
 #endif
 #ifdef WARPX_DIM_RZ
-        m_structs[i].pos(0) = x;
+        m_x[i] = x;
         m_theta[i] = y;
-        m_structs[i].pos(1) = z;
+        m_z[i] = z;
 #elif WARPX_DIM_3D
-        m_structs[i].pos(0) = x;
-        m_structs[i].pos(1) = y;
-        m_structs[i].pos(2) = z;
+        m_x[i] = x;
+        m_y[i] = y;
+        m_z[i] = z;
 #elif WARPX_DIM_XZ
-        m_structs[i].pos(0) = x;
-        m_structs[i].pos(1) = z;
+        m_x[i] = x;
+        m_z[i] = z;
 #else
-        m_structs[i].pos(0) = z;
+        m_z[i] = z;
 #endif
     }
 };
diff --git a/Source/Particles/Resampling/LevelingThinning.cpp b/Source/Particles/Resampling/LevelingThinning.cpp
index 680e33ebe6a..40f75c76eab 100644
--- a/Source/Particles/Resampling/LevelingThinning.cpp
+++ b/Source/Particles/Resampling/LevelingThinning.cpp
@@ -60,8 +60,7 @@ void LevelingThinning::operator() (WarpXParIter& pti, const int lev,
     auto& ptile = pc->ParticlesAt(lev, pti);
     auto& soa = ptile.GetStructOfArrays();
     amrex::ParticleReal * const AMREX_RESTRICT w = soa.GetRealData(PIdx::w).data();
-    WarpXParticleContainer::ParticleType * const AMREX_RESTRICT
-                                 particle_ptr = ptile.GetArrayOfStructs()().data();
+    auto * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data();
 
     // Using this function means that we must loop over the cells in the ParallelFor. In the case
     // of the leveling thinning algorithm, it would have possibly been more natural and more
@@ -114,7 +113,7 @@ void LevelingThinning::operator() (WarpXParIter& pti, const int lev,
                 // Remove particle with probability 1 - particle_weight/level_weight
                 if (random_number > w[indices[i]]/level_weight)
                 {
-                    particle_ptr[indices[i]].id() = -1;
+                    amrex::ParticleIDWrapper{idcpu[indices[i]]} = -1;
                 }
                 // Set particle weight to level weight otherwise
                 else
diff --git a/Source/Particles/Sorting/Partition.cpp b/Source/Particles/Sorting/Partition.cpp
index 58511cfd5e7..58e3450f47d 100644
--- a/Source/Particles/Sorting/Partition.cpp
+++ b/Source/Particles/Sorting/Partition.cpp
@@ -61,7 +61,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers(
     // Initialize temporary arrays
     Gpu::DeviceVector<int> inexflag;
     inexflag.resize(np);
-    Gpu::DeviceVector<long> pid;
+    Gpu::DeviceVector<int> pid;
     pid.resize(np);
 
     // First, partition particles into the larger buffer
@@ -109,7 +109,7 @@ PhysicalParticleContainer::PartitionParticlesInBuffers(
             // - For each particle in the large buffer, find whether it is in
             // the smaller buffer, by looking up the mask. Store the answer in `inexflag`.
             amrex::ParallelFor( np - n_fine,
-               fillBufferFlagRemainingParticles(pti, bmasks, inexflag, Geom(lev), pid, n_fine) );
+               fillBufferFlagRemainingParticles(pti, bmasks, inexflag, Geom(lev), pid, int(n_fine)) );
             auto *const sep2 = stablePartition( sep, pid.end(), inexflag );
 
             if (bmasks == gather_masks) {
diff --git a/Source/Particles/Sorting/SortingUtils.H b/Source/Particles/Sorting/SortingUtils.H
index ac2c63e88f8..ba7761bf48a 100644
--- a/Source/Particles/Sorting/SortingUtils.H
+++ b/Source/Particles/Sorting/SortingUtils.H
@@ -12,6 +12,7 @@
 
 #include <AMReX_Gpu.H>
 #include <AMReX_Partition.H>
+#include <AMReX_ParticleUtil.H>
 
 
 /** \brief Fill the elements of the input vector with consecutive integer,
@@ -19,7 +20,7 @@
  *
  * \param[inout] v Vector of integers, to be filled by this routine
  */
-void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector<long>& v );
+void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector<int>& v );
 
 /** \brief Find the indices that would reorder the elements of `predicate`
  * so that the elements with non-zero value precede the other elements
@@ -41,7 +42,7 @@ ForwardIterator stablePartition(ForwardIterator const index_begin,
     int const* AMREX_RESTRICT predicate_ptr = predicate.dataPtr();
     int N = static_cast<int>(std::distance(index_begin, index_end));
     auto num_true = amrex::StablePartition(&(*index_begin), N,
-        [predicate_ptr] AMREX_GPU_DEVICE (long i) { return predicate_ptr[i]; });
+        [predicate_ptr] AMREX_GPU_DEVICE (int i) { return predicate_ptr[i]; });
 
     ForwardIterator sep = index_begin;
     std::advance(sep, num_true);
@@ -49,7 +50,7 @@ ForwardIterator stablePartition(ForwardIterator const index_begin,
     // On CPU: Use std library
     ForwardIterator const sep = std::stable_partition(
         index_begin, index_end,
-        [&predicate](long i) { return predicate[i]; }
+        [&predicate](int i) { return predicate[i]; }
     );
 #endif
     return sep;
@@ -88,7 +89,7 @@ class fillBufferFlag
             // Extract simple structure that can be used directly on the GPU
             m_domain{geom.Domain()},
             m_inexflag_ptr{inexflag.dataPtr()},
-            m_particles{pti.GetArrayOfStructs().data()},
+            m_ptd{pti.GetParticleTile().getConstParticleTileData()},
             m_buffer_mask{(*bmasks)[pti].array()}
         {
             for (int idim=0; idim<AMREX_SPACEDIM; idim++) {
@@ -99,11 +100,9 @@ class fillBufferFlag
 
 
         AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-        void operator()( const long i ) const {
-            // Select a particle
-            auto const& p = m_particles[i];
+        void operator()( const int i ) const {
             // Find the index of the cell where this particle is located
-            amrex::IntVect const iv = amrex::getParticleCell( p,
+            amrex::IntVect const iv = amrex::getParticleCell( m_ptd, i,
                                 m_prob_lo, m_inv_cell_size, m_domain );
             // Find the value of the buffer flag in this cell and
             // store it at the corresponding particle position in the array `inexflag`
@@ -113,7 +112,7 @@ class fillBufferFlag
     private:
         amrex::Box m_domain;
         int* m_inexflag_ptr;
-        WarpXParticleContainer::ParticleType const* m_particles;
+        WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType m_ptd;
         amrex::Array4<int const> m_buffer_mask;
         amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> m_prob_lo;
         amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> m_inv_cell_size;
@@ -141,12 +140,12 @@ class fillBufferFlagRemainingParticles
                         amrex::iMultiFab const* bmasks,
                         amrex::Gpu::DeviceVector<int>& inexflag,
                         amrex::Geometry const& geom,
-                        amrex::Gpu::DeviceVector<long> const& particle_indices,
-                        long const start_index ) :
+                        amrex::Gpu::DeviceVector<int> const& particle_indices,
+                        int start_index ) :
             m_domain{geom.Domain()},
             // Extract simple structure that can be used directly on the GPU
             m_inexflag_ptr{inexflag.dataPtr()},
-            m_particles{pti.GetArrayOfStructs().data()},
+            m_ptd{pti.GetParticleTile().getConstParticleTileData()},
             m_buffer_mask{(*bmasks)[pti].array()},
             m_start_index{start_index},
             m_indices_ptr{particle_indices.dataPtr()}
@@ -159,11 +158,11 @@ class fillBufferFlagRemainingParticles
 
 
         AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-        void operator()( const long i ) const {
+        void operator()( const int i ) const {
             // Select a particle
-            auto const& p = m_particles[m_indices_ptr[i+m_start_index]];
+            auto const j = m_indices_ptr[i+m_start_index];
             // Find the index of the cell where this particle is located
-            amrex::IntVect const iv = amrex::getParticleCell( p,
+            amrex::IntVect const iv = amrex::getParticleCell( m_ptd, j,
                                 m_prob_lo, m_inv_cell_size, m_domain );
             // Find the value of the buffer flag in this cell and
             // store it at the corresponding particle position in the array `inexflag`
@@ -175,10 +174,10 @@ class fillBufferFlagRemainingParticles
         amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> m_inv_cell_size;
         amrex::Box m_domain;
         int* m_inexflag_ptr;
-        WarpXParticleContainer::ParticleType const* m_particles;
+        WarpXParticleContainer::ParticleTileType::ConstParticleTileDataType const m_ptd;
         amrex::Array4<int const> m_buffer_mask;
-        long const m_start_index;
-        long const* m_indices_ptr;
+        int const m_start_index;
+        int const* m_indices_ptr;
 };
 
 /** \brief Functor that copies the elements of `src` into `dst`,
@@ -195,7 +194,7 @@ class copyAndReorder
         copyAndReorder(
             amrex::Gpu::DeviceVector<T> const& src,
             amrex::Gpu::DeviceVector<T>& dst,
-            amrex::Gpu::DeviceVector<long> const& indices ):
+            amrex::Gpu::DeviceVector<int> const& indices ):
                 // Extract simple structure that can be used directly on the GPU
                 m_src_ptr{src.dataPtr()},
                 m_dst_ptr{dst.dataPtr()},
@@ -203,14 +202,14 @@ class copyAndReorder
         {}
 
         AMREX_GPU_DEVICE AMREX_FORCE_INLINE
-        void operator()( const long ip ) const {
+        void operator()( const int ip ) const {
             m_dst_ptr[ip] = m_src_ptr[ m_indices_ptr[ip] ];
         }
 
     private:
         T const* m_src_ptr;
         T* m_dst_ptr;
-        long const* m_indices_ptr;
+        int const* m_indices_ptr;
 };
 
 #endif // WARPX_PARTICLES_SORTING_SORTINGUTILS_H_
diff --git a/Source/Particles/Sorting/SortingUtils.cpp b/Source/Particles/Sorting/SortingUtils.cpp
index 699119e8e18..cd4b6a13c76 100644
--- a/Source/Particles/Sorting/SortingUtils.cpp
+++ b/Source/Particles/Sorting/SortingUtils.cpp
@@ -8,7 +8,7 @@
 
 #include "SortingUtils.H"
 
-void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector<long>& v )
+void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector<int>& v )
 {
 #ifdef AMREX_USE_GPU
     // On GPU: Use amrex
diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H
index 33aa71d1c7d..7d2d5619da9 100644
--- a/Source/Particles/WarpXParticleContainer.H
+++ b/Source/Particles/WarpXParticleContainer.H
@@ -49,10 +49,10 @@
 
 
 class WarpXParIter
-    : public amrex::ParIter<0,0,PIdx::nattribs>
+    : public amrex::ParIterSoA<PIdx::nattribs, 0>
 {
 public:
-    using amrex::ParIter<0,0,PIdx::nattribs>::ParIter;
+    using amrex::ParIterSoA<PIdx::nattribs, 0>::ParIterSoA;
 
     WarpXParIter (ContainerType& pc, int level);
 
@@ -89,13 +89,14 @@ public:
  * particle container classes (that store a collection of particles) derive. Derived
  * classes can be used for plasma particles, photon particles, or non-physical
  * particles (e.g., for the laser antenna).
- * It derives from amrex::ParticleContainer<0,0,PIdx::nattribs>, where the
- * template arguments stand for the number of int and amrex::Real SoA and AoS
- * data in amrex::Particle.
- *  - AoS amrex::Real: x, y, z (default), 0 additional (first template
- *    parameter)
- *  - AoS int: id, cpu (default), 0 additional (second template parameter)
- *  - SoA amrex::Real: PIdx::nattribs (third template parameter), see PIdx for
+ * It derives from amrex::ParticleContainerPureSoA<PIdx::nattribs>, where the
+ * template arguments stand for the number of int and amrex::Real SoA
+ * data in amrex::SoAParticle.
+ *  - SoA amrex::Real: positions x, y, z, momentum ux, uy, uz, ... see PIdx for details;
+ *                     more can be added at runtime
+ *  - SoA int: 0 attributes by default, but can be added at runtime
+ *  - SoA uint64_t: idcpu, a global 64bit index, with a 40bit local id and a 24bit cpu id
+ *    (both set at creation)
  * the list.
  *
  * WarpXParticleContainer contains the main functions for initialization,
@@ -164,11 +165,9 @@ public:
      * class.
      */
     virtual void DefaultInitializeRuntimeAttributes (
-                        amrex::ParticleTile<amrex::Particle<NStructReal, NStructInt>,
-                                            NArrayReal, NArrayInt,
-                                            amrex::PinnedArenaAllocator>& pinned_tile,
-                        int n_external_attr_real,
-                        int n_external_attr_int) = 0;
+        typename ContainerLike<amrex::PinnedArenaAllocator>::ParticleTileType& pinned_tile,
+        int n_external_attr_real,
+        int n_external_attr_int) = 0;
 
     ///
     /// This pushes the particle positions by one half time step.
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index a395198e361..eb590f893e9 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -75,13 +75,13 @@
 using namespace amrex;
 
 WarpXParIter::WarpXParIter (ContainerType& pc, int level)
-    : amrex::ParIter<0,0,PIdx::nattribs>(pc, level,
+    : amrex::ParIterSoA<PIdx::nattribs, 0>(pc, level,
              MFItInfo().SetDynamic(WarpX::do_dynamic_scheduling))
 {
 }
 
 WarpXParIter::WarpXParIter (ContainerType& pc, int level, MFItInfo& info)
-    : amrex::ParIter<0,0,PIdx::nattribs>(pc, level,
+    : amrex::ParIterSoA<PIdx::nattribs, 0>(pc, level,
                    info.SetDynamic(WarpX::do_dynamic_scheduling))
 {
 }
@@ -198,52 +198,54 @@ WarpXParticleContainer::AddNParticles (int /*lev*/, long n,
     // Redistribute() will move them to proper places.
     auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0);
 
-    using PinnedTile = amrex::ParticleTile<Particle<NStructReal, NStructInt>,
-                                           NArrayReal, NArrayInt,
-                                           amrex::PinnedArenaAllocator>;
+    using PinnedTile = typename ContainerLike<amrex::PinnedArenaAllocator>::ParticleTileType;
     PinnedTile pinned_tile;
     pinned_tile.define(NumRuntimeRealComps(), NumRuntimeIntComps());
 
     const std::size_t np = iend-ibegin;
 
 #ifdef WARPX_DIM_RZ
+    amrex::Vector<amrex::ParticleReal> r(np);
     amrex::Vector<amrex::ParticleReal> theta(np);
 #endif
 
     for (auto i = ibegin; i < iend; ++i)
     {
-        ParticleType p;
-        if (id==-1)
-        {
-            p.id() = ParticleType::NextID();
+        auto & idcpu_data = pinned_tile.GetStructOfArrays().GetIdCPUData();
+        idcpu_data.push_back(0);
+        if (id==-1) {
+            amrex::ParticleIDWrapper{idcpu_data.back()} = ParticleType::NextID();
         } else {
-            p.id() = id;
+            amrex::ParticleIDWrapper{idcpu_data.back()} = id;
         }
-        p.cpu() = amrex::ParallelDescriptor::MyProc();
+        amrex::ParticleCPUWrapper(idcpu_data.back()) = ParallelDescriptor::MyProc();
+
+#ifdef WARPX_DIM_RZ
+        r[i-ibegin] = std::sqrt(x[i]*x[i] + y[i]*y[i]);
+        theta[i-ibegin] = std::atan2(y[i], x[i]);
+#endif
+    }
+
+    if (np > 0)
+    {
 #if defined(WARPX_DIM_3D)
-        p.pos(0) = x[i];
-        p.pos(1) = y[i];
-        p.pos(2) = z[i];
+        pinned_tile.push_back_real(PIdx::x, x.data() + ibegin, x.data() + iend);
+        pinned_tile.push_back_real(PIdx::y, y.data() + ibegin, y.data() + iend);
+        pinned_tile.push_back_real(PIdx::z, z.data() + ibegin, z.data() + iend);
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
         amrex::ignore_unused(y);
 #ifdef WARPX_DIM_RZ
-        theta[i-ibegin] = std::atan2(y[i], x[i]);
-        p.pos(0) = std::sqrt(x[i]*x[i] + y[i]*y[i]);
+        pinned_tile.push_back_real(PIdx::x, r.data(), r.data() + np);
 #else
-        p.pos(0) = x[i];
+        pinned_tile.push_back_real(PIdx::x, x.data() + ibegin, x.data() + iend);
 #endif
-        p.pos(1) = z[i];
+        pinned_tile.push_back_real(PIdx::z, z.data() + ibegin, z.data() + iend);
 #else //AMREX_SPACEDIM == 1
         amrex::ignore_unused(x,y);
-        p.pos(0) = z[i];
+        pinned_tile.push_back_real(PIdx::z, z.data() + ibegin, z.data() + iend);
 #endif
 
-        pinned_tile.push_back(p);
-    }
-
-    if (np > 0)
-    {
-        pinned_tile.push_back_real(PIdx::w , attr_real[0].data() + ibegin, attr_real[0].data() + iend);
+        pinned_tile.push_back_real(PIdx::w, attr_real[0].data() + ibegin, attr_real[0].data() + iend);
         pinned_tile.push_back_real(PIdx::ux, ux.data() + ibegin, ux.data() + iend);
         pinned_tile.push_back_real(PIdx::uy, uy.data() + ibegin, uy.data() + iend);
         pinned_tile.push_back_real(PIdx::uz, uz.data() + ibegin, uz.data() + iend);
@@ -476,15 +478,14 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti,
 
         //sort particles by bin
         WARPX_PROFILE_VAR_START(blp_sort);
-        amrex::DenseBins<ParticleType> bins;
+        amrex::DenseBins<ParticleTileType::ParticleTileDataType> bins;
         {
             auto& ptile = ParticlesAt(lev, pti);
-            auto& aos = ptile.GetArrayOfStructs();
-            auto *pstruct_ptr = aos().dataPtr();
+            auto ptd = ptile.getParticleTileData();
 
             const int ntiles = numTilesInBox(box, true, bin_size);
 
-            bins.build(ptile.numParticles(), pstruct_ptr, ntiles,
+            bins.build(ptile.numParticles(), ptd, ntiles,
                     [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) -> unsigned int
                     {
                         Box tbox;
@@ -947,7 +948,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp,
 
         // HACK - sort particles by bin here.
         WARPX_PROFILE_VAR_START(blp_sort);
-        amrex::DenseBins<ParticleType> bins;
+        amrex::DenseBins<ParticleTileType::ParticleTileDataType> bins;
         {
             const Geometry& geom = Geom(lev);
             const auto dxi = geom.InvCellSizeArray();
@@ -955,16 +956,15 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp,
             const auto domain = geom.Domain();
 
             auto& ptile = ParticlesAt(lev, pti);
-            auto& aos   = ptile.GetArrayOfStructs();
-            auto *pstruct_ptr = aos().dataPtr();
+            auto ptd = ptile.getParticleTileData();
 
             Box box = pti.validbox();
             box.grow(ng_rho);
             const amrex::IntVect bin_size = WarpX::shared_tilesize;
             const int ntiles = numTilesInBox(box, true, bin_size);
 
-            bins.build(ptile.numParticles(), pstruct_ptr, ntiles,
-                       [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) -> unsigned int
+            bins.build(ptile.numParticles(), ptd, ntiles,
+                       [=] AMREX_GPU_HOST_DEVICE (ParticleType const & p) -> unsigned int
                        {
                            Box tbx;
                            auto iv = getParticleCell(p, plo, dxi, domain);
@@ -984,8 +984,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp,
             const auto domain = geom.Domain();
 
             auto& ptile = ParticlesAt(lev, pti);
-            auto& aos   = ptile.GetArrayOfStructs();
-            auto *pstruct_ptr = aos().dataPtr();
+            auto ptd = ptile.getParticleTileData();
 
             Box box = pti.validbox();
             box.grow(ng_rho);
@@ -999,9 +998,10 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp,
                                    const auto bin_start = offsets_ptr[ibin];
                                    const auto bin_stop = offsets_ptr[ibin+1];
                                    if (bin_start < bin_stop) {
-                                       auto p = pstruct_ptr[permutation[bin_start]];
+                                       // static_cast until https://github.com/AMReX-Codes/amrex/pull/3684
+                                       auto const i = static_cast<int>(permutation[bin_start]);
                                        Box tbx;
-                                       auto iv = getParticleCell(p, plo, dxi, domain);
+                                       auto iv = getParticleCell(ptd, i, plo, dxi, domain);
                                        AMREX_ASSERT(box.contains(iv));
                                        [[maybe_unused]] auto tid = getTileIndex(iv, box, true, bin_size, tbx);
                                        AMREX_ASSERT(tid == ibin);
@@ -1490,10 +1490,10 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p,
     // Tag particle if goes to higher level.
     // It will be split later in the loop
     if (pld.m_lev == lev+1
-        and p.id() != NoSplitParticleID
+        and p.id() != amrex::LongParticleIds::NoSplitParticleID
         and p.id() >= 0)
     {
-        p.id() = DoSplitParticleID;
+        p.id() = amrex::LongParticleIds::DoSplitParticleID;
     }
 
     if (pld.m_lev == lev-1){
@@ -1532,9 +1532,9 @@ WarpXParticleContainer::ApplyBoundaryConditions (){
             const Real zmax = Geom(lev).ProbHi(WARPX_ZINDEX);
 
             ParticleTileType& ptile = ParticlesAt(lev, pti);
-            ParticleType * const pp = ptile.GetArrayOfStructs()().data();
 
             auto& soa = ptile.GetStructOfArrays();
+            uint64_t * const AMREX_RESTRICT idcpu = soa.GetIdCPUData().data();
             amrex::ParticleReal * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data();
             amrex::ParticleReal * const AMREX_RESTRICT uy = soa.GetRealData(PIdx::uy).data();
             amrex::ParticleReal * const AMREX_RESTRICT uz = soa.GetRealData(PIdx::uz).data();
@@ -1543,10 +1543,9 @@ WarpXParticleContainer::ApplyBoundaryConditions (){
             amrex::ParallelForRNG(
                 pti.numParticles(),
                 [=] AMREX_GPU_DEVICE (long i, amrex::RandomEngine const& engine) {
-                    ParticleType& p = pp[i];
-
                     // skip particles that are already flagged for removal
-                    if (p.id() < 0) { return; }
+                    auto const id = amrex::ParticleIDWrapper{idcpu[i]};
+                    if (id < 0) { return; }
 
                     ParticleReal x, y, z;
                     GetPosition.AsStored(i, x, y, z);
@@ -1568,7 +1567,7 @@ WarpXParticleContainer::ApplyBoundaryConditions (){
                                                               boundary_conditions, engine);
 
                     if (particle_lost) {
-                        p.id() = -p.id();
+                        amrex::ParticleIDWrapper{idcpu[i]} = -id;
                     } else {
                         SetPosition.AsStored(i, x, y, z);
                     }
diff --git a/Source/Python/Particles/ParticleBoundaryBuffer.cpp b/Source/Python/Particles/ParticleBoundaryBuffer.cpp
index 2a35faece9b..b04ac75e600 100644
--- a/Source/Python/Particles/ParticleBoundaryBuffer.cpp
+++ b/Source/Python/Particles/ParticleBoundaryBuffer.cpp
@@ -10,13 +10,13 @@
 
 namespace warpx {
     class BoundaryBufferParIter
-        : public amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>
+        : public amrex::ParIterSoA<PIdx::nattribs, 0, amrex::PinnedArenaAllocator>
     {
     public:
-        using amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>::ParIter;
+        using amrex::ParIterSoA<PIdx::nattribs, 0, amrex::PinnedArenaAllocator>::ParIterSoA;
 
         BoundaryBufferParIter(ContainerType& pc, int level) :
-            amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>(pc, level) {}
+            amrex::ParIterSoA<PIdx::nattribs, 0, amrex::PinnedArenaAllocator>(pc, level) {}
     };
 }
 
@@ -24,9 +24,9 @@ void init_BoundaryBufferParIter (py::module& m)
 {
     py::class_<
         warpx::BoundaryBufferParIter,
-        amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>
+        amrex::ParIterSoA<PIdx::nattribs, 0, amrex::PinnedArenaAllocator>
     >(m, "BoundaryBufferParIter")
-        .def(py::init<amrex::ParIter<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>::ContainerType&, int>(),
+        .def(py::init<amrex::ParIterSoA<PIdx::nattribs, 0, amrex::PinnedArenaAllocator>::ContainerType&, int>(),
             py::arg("particle_container"), py::arg("level")
         )
     ;
diff --git a/Source/Python/Particles/PinnedMemoryParticleContainer.cpp b/Source/Python/Particles/PinnedMemoryParticleContainer.cpp
index 600d56a62c9..d4f6a422dbe 100644
--- a/Source/Python/Particles/PinnedMemoryParticleContainer.cpp
+++ b/Source/Python/Particles/PinnedMemoryParticleContainer.cpp
@@ -13,6 +13,6 @@ void init_PinnedMemoryParticleContainer (py::module& m)
 {
     py::class_<
         PinnedMemoryParticleContainer,
-        amrex::ParticleContainer<0,0,PIdx::nattribs,0,amrex::PinnedArenaAllocator>
+        amrex::ParticleContainerPureSoA<PIdx::nattribs, 0, amrex::PinnedArenaAllocator>
     > pmpc (m, "PinnedMemoryParticleContainer");
 }
diff --git a/Source/Python/Particles/WarpXParticleContainer.cpp b/Source/Python/Particles/WarpXParticleContainer.cpp
index 1473a750941..07793a373f3 100644
--- a/Source/Python/Particles/WarpXParticleContainer.cpp
+++ b/Source/Python/Particles/WarpXParticleContainer.cpp
@@ -12,11 +12,11 @@
 void init_WarpXParIter (py::module& m)
 {
     py::class_<
-        WarpXParIter, amrex::ParIter<0,0,PIdx::nattribs>
+        WarpXParIter, amrex::ParIterSoA<PIdx::nattribs, 0>
     >(m, "WarpXParIter")
-        .def(py::init<amrex::ParIter<0,0,PIdx::nattribs>::ContainerType&, int>(),
+        .def(py::init<amrex::ParIterSoA<PIdx::nattribs, 0>::ContainerType&, int>(),
             py::arg("particle_container"), py::arg("level"))
-        .def(py::init<amrex::ParIter<0,0,PIdx::nattribs>::ContainerType&, int, amrex::MFItInfo&>(),
+        .def(py::init<amrex::ParIterSoA<PIdx::nattribs, 0>::ContainerType&, int, amrex::MFItInfo&>(),
             py::arg("particle_container"), py::arg("level"),
             py::arg("info"))
     ;
@@ -26,11 +26,11 @@ void init_WarpXParticleContainer (py::module& m)
 {
     py::class_<
         WarpXParticleContainer,
-        amrex::ParticleContainer<0, 0, PIdx::nattribs, 0>
+        amrex::ParticleContainerPureSoA<PIdx::nattribs, 0>
     > wpc (m, "WarpXParticleContainer");
     wpc
         .def("add_real_comp",
-            [](WarpXParticleContainer& pc, const std::string& name, bool const comm) { pc.AddRealComp(name, comm); },
+            [](WarpXParticleContainer& pc, const std::string& name, bool comm) { pc.AddRealComp(name, comm); },
             py::arg("name"), py::arg("comm")
         )
         .def("add_n_particles",
@@ -93,6 +93,14 @@ void init_WarpXParticleContainer (py::module& m)
             },
             py::arg("comp_name")
         )
+        .def("get_icomp_index",
+            [](WarpXParticleContainer& pc, std::string comp_name)
+            {
+                auto particle_comps = pc.getParticleiComps();
+                return particle_comps.at(comp_name);
+            },
+            py::arg("comp_name")
+        )
         .def("num_local_tiles_at_level",
             &WarpXParticleContainer::numLocalTilesAtLevel,
             py::arg("level")
diff --git a/Source/Utils/ParticleUtils.H b/Source/Utils/ParticleUtils.H
index b04176d4d83..7e3c89228ea 100644
--- a/Source/Utils/ParticleUtils.H
+++ b/Source/Utils/ParticleUtils.H
@@ -28,9 +28,10 @@ namespace ParticleUtils {
      * @param[in] mfi the MultiFAB iterator.
      * @param[in] ptile the particle tile.
      */
-    amrex::DenseBins<WarpXParticleContainer::ParticleType>
-    findParticlesInEachCell(int lev, amrex::MFIter const& mfi,
-                             WarpXParticleContainer::ParticleTileType const& ptile);
+    amrex::DenseBins<typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType>
+    findParticlesInEachCell (int lev,
+                             amrex::MFIter const & mfi,
+                             WarpXParticleContainer::ParticleTileType & ptile);
 
     /**
      * \brief Return (relativistic) particle energy given velocity and mass.
diff --git a/Source/Utils/ParticleUtils.cpp b/Source/Utils/ParticleUtils.cpp
index 60e04f12b86..b8207b61fa0 100644
--- a/Source/Utils/ParticleUtils.cpp
+++ b/Source/Utils/ParticleUtils.cpp
@@ -22,24 +22,28 @@
 #include <AMReX_REAL.H>
 #include <AMReX_SPACE.H>
 
-namespace ParticleUtils {
+namespace ParticleUtils
+{
 
     using namespace amrex;
+
     // Define shortcuts for frequently-used type names
-    using ParticleType = WarpXParticleContainer::ParticleType;
-    using ParticleTileType = WarpXParticleContainer::ParticleTileType;
-    using ParticleBins = DenseBins<ParticleType>;
-    using index_type = ParticleBins::index_type;
+    using ParticleType = typename WarpXParticleContainer::ParticleType;
+    using ParticleTileType = typename WarpXParticleContainer::ParticleTileType;
+    using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType;
+    using ParticleBins = DenseBins<ParticleTileDataType>;
+    using index_type = typename ParticleBins::index_type;
 
     /* Find the particles and count the particles that are in each cell.
        Note that this does *not* rearrange particle arrays */
     ParticleBins
-    findParticlesInEachCell( int const lev, MFIter const& mfi,
-                             ParticleTileType const& ptile) {
+    findParticlesInEachCell (int lev,
+                             MFIter const & mfi,
+                             ParticleTileType & ptile) {
 
         // Extract particle structures for this tile
         int const np = ptile.numParticles();
-        ParticleType const* particle_ptr = ptile.GetArrayOfStructs()().data();
+        auto ptd = ptile.getParticleTileData();
 
         // Extract box properties
         Geometry const& geom = WarpX::GetInstance().Geom(lev);
@@ -51,9 +55,9 @@ namespace ParticleUtils {
         // Find particles that are in each cell;
         // results are stored in the object `bins`.
         ParticleBins bins;
-        bins.build(np, particle_ptr, cbx,
+        bins.build(np, ptd, cbx,
             // Pass lambda function that returns the cell index
-            [=] AMREX_GPU_DEVICE (const ParticleType& p) noexcept
+            [=] AMREX_GPU_DEVICE (ParticleType const & p) noexcept -> amrex::IntVect
             {
                 return IntVect{AMREX_D_DECL(
                                    static_cast<int>((p.pos(0)-plo[0])*dxi[0] - lo.x),
@@ -64,4 +68,4 @@ namespace ParticleUtils {
         return bins;
     }
 
-}
+} // namespace ParticleUtils
diff --git a/Source/ablastr/particles/IndexHandling.H b/Source/ablastr/particles/IndexHandling.H
deleted file mode 100644
index 0ad5ca60446..00000000000
--- a/Source/ablastr/particles/IndexHandling.H
+++ /dev/null
@@ -1,41 +0,0 @@
-/* Copyright 2019-2022 Axel Huebl
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
-#ifndef ABLASTR_INDEX_HANDLING_H
-#define ABLASTR_INDEX_HANDLING_H
-
-#include <cstdint>
-
-
-namespace ablastr::particles {
-
-    /** A helper function to derive a globally unique particle ID
-     *
-     * @param[in] id  AMReX particle ID (on local cpu/rank), AoS .id
-     * @param[in] cpu AMReX particle CPU (rank) at creation of the particle, AoS .cpu
-     * @return global particle ID that is unique and permanent in the whole simulation
-     */
-    constexpr uint64_t
-    localIDtoGlobal (int const id, int const cpu)
-    {
-        static_assert(sizeof(int) * 2u <= sizeof(uint64_t),
-                      "int size might cause collisions in global IDs");
-        // implementation:
-        //   - we cast both 32-bit (or smaller) ints to a 64bit unsigned int
-        //   - this will leave half of the "upper range" bits in the 64bit unsigned int zeroed out
-        //     because the corresponding (extended) value range was not part of the value range in
-        //     the int representation
-        //   - we bit-shift the cpu into the upper half of zero bits in the 64 bit unsigned int
-        //     (imagine this step as "adding a per-cpu/rank offset to the local integers")
-        //   - then we add this offset
-        //     note: the add is expressed as bitwise OR (|) since this saves us from writing
-        //           brackets for operator precedence between + and <<
-        return uint64_t(id) | uint64_t(cpu) << 32u;
-    }
-
-} // namespace ablastr::particles
-
-#endif // ABLASTR_INDEX_HANDLING_H
diff --git a/Source/ablastr/particles/ParticleMoments.H b/Source/ablastr/particles/ParticleMoments.H
index e45fb574cce..b648ccb28aa 100644
--- a/Source/ablastr/particles/ParticleMoments.H
+++ b/Source/ablastr/particles/ParticleMoments.H
@@ -35,7 +35,7 @@ namespace particles {
             amrex::ParticleReal, amrex::ParticleReal>
     MinAndMaxPositions (T_PC const & pc)
     {
-        using PType = typename T_PC::SuperParticleType;
+        using ConstParticleTileDataType = typename T_PC::ParticleTileType::ConstParticleTileDataType;
 
         // Get min and max for the local rank
         amrex::ReduceOps<
@@ -46,11 +46,11 @@ namespace particles {
                               amrex::ParticleReal, amrex::ParticleReal, amrex::ParticleReal>
         >(
             pc,
-            [=] AMREX_GPU_DEVICE(PType const & p) noexcept
+            [=] AMREX_GPU_DEVICE(const ConstParticleTileDataType& ptd, const int i) noexcept
             {
-                amrex::ParticleReal const x = p.pos(0);
-                amrex::ParticleReal const y = p.pos(1);
-                amrex::ParticleReal const z = p.pos(2);
+                const amrex::ParticleReal x = ptd.rdata(0)[i];
+                const amrex::ParticleReal y = ptd.rdata(1)[i];
+                const amrex::ParticleReal z = ptd.rdata(2)[i];
 
                 return amrex::makeTuple(x, y, z, x, y, z);
             },
@@ -90,7 +90,8 @@ namespace particles {
             amrex::ParticleReal, amrex::ParticleReal>
     MeanAndStdPositions (T_PC const & pc)
     {
-        using PType = typename T_PC::SuperParticleType;
+
+        using ConstParticleTileDataType = typename T_PC::ParticleTileType::ConstParticleTileDataType;
 
         amrex::ReduceOps<
             amrex::ReduceOpSum, amrex::ReduceOpSum, amrex::ReduceOpSum,
@@ -103,12 +104,14 @@ namespace particles {
                 amrex::ParticleReal>
         >(
             pc,
-            [=] AMREX_GPU_DEVICE(const PType& p) noexcept
+            [=] AMREX_GPU_DEVICE(const ConstParticleTileDataType& ptd, const int i) noexcept
             {
-                amrex::ParticleReal const x = p.pos(0);
-                amrex::ParticleReal const y = p.pos(1);
-                amrex::ParticleReal const z = p.pos(2);
-                amrex::ParticleReal const w = p.rdata(T_RealSoAWeight);
+
+                const amrex::ParticleReal x = ptd.rdata(0)[i];
+                const amrex::ParticleReal y = ptd.rdata(1)[i];
+                const amrex::ParticleReal z = ptd.rdata(2)[i];
+
+                const amrex::ParticleReal w = ptd.rdata(T_RealSoAWeight)[i];
 
                 return amrex::makeTuple(x, x*x, y, y*y, z, z*z, w);
             },