diff --git a/.github/workflows/dependencies/hip.sh b/.github/workflows/dependencies/hip.sh
index 2225e670bb0..ae897786acd 100755
--- a/.github/workflows/dependencies/hip.sh
+++ b/.github/workflows/dependencies/hip.sh
@@ -12,12 +12,24 @@ set -eu -o pipefail
 #   failed files the given number of times.
 echo 'Acquire::Retries "3";' | sudo tee /etc/apt/apt.conf.d/80-retries
 
-# Ref.: https://rocmdocs.amd.com/en/latest/Installation_Guide/Installation-Guide.html#ubuntu
-wget -q -O - https://repo.radeon.com/rocm/rocm.gpg.key \
-  | sudo apt-key add -
-echo 'deb [arch=amd64] https://repo.radeon.com/rocm/apt/debian/ ubuntu main' \
-  | sudo tee /etc/apt/sources.list.d/rocm.list
+# Ref.: https://rocm.docs.amd.com/projects/install-on-linux/en/latest/how-to/native-install/ubuntu.html
+
+# Make the directory if it doesn't exist yet.
+# This location is recommended by the distribution maintainers.
+sudo mkdir --parents --mode=0755 /etc/apt/keyrings
+
+# Download the key, convert the signing-key to a full
+# keyring required by apt and store in the keyring directory
+wget https://repo.radeon.com/rocm/rocm.gpg.key -O - | \
+    gpg --dearmor | sudo tee /etc/apt/keyrings/rocm.gpg > /dev/null
 
+curl -O https://repo.radeon.com/rocm/rocm.gpg.key
+sudo apt-key add rocm.gpg.key
+
+source /etc/os-release # set UBUNTU_CODENAME: focal or jammy or ...
+
+echo "deb [arch=amd64] https://repo.radeon.com/rocm/apt/${1-latest} ${UBUNTU_CODENAME} main" \
+  | sudo tee /etc/apt/sources.list.d/rocm.list
 echo 'export PATH=/opt/rocm/llvm/bin:/opt/rocm/bin:/opt/rocm/profiler/bin:/opt/rocm/opencl/bin:$PATH' \
   | sudo tee -a /etc/profile.d/rocm.sh
 
diff --git a/Docs/source/highlights.rst b/Docs/source/highlights.rst
index 91c99ee50ca..b243f62cd97 100644
--- a/Docs/source/highlights.rst
+++ b/Docs/source/highlights.rst
@@ -25,9 +25,9 @@ Scientific works in laser-plasma and beam-plasma acceleration.
    `DOI:10.1103/PhysRevResearch.5.033112 <https://doi.org/10.1103/PhysRevResearch.5.033112>`__
 
 #. Sandberg R T, Lehe R, Mitchell C E, Garten M, Myers A, Qiang J, Vay J-L and Huebl A.
-  **Synthesizing Particle-in-Cell Simulations Through Learning and GPU Computing for Hybrid Particle Accelerator Beamlines**.
-  Proc. of Platform for Advanced Scientific Computing (PASC'24), *submitted*, 2024.
-  `preprint <http://arxiv.org/abs/2402.17248>__`
+   **Synthesizing Particle-in-Cell Simulations Through Learning and GPU Computing for Hybrid Particle Accelerator Beamlines**.
+   Proc. of Platform for Advanced Scientific Computing (PASC'24), *submitted*, 2024.
+   `preprint <http://arxiv.org/abs/2402.17248>`__
 
 #. Sandberg R T, Lehe R, Mitchell C E, Garten M, Qiang J, Vay J-L and Huebl A.
    **Hybrid Beamline Element ML-Training for Surrogates in the ImpactX Beam-Dynamics Code**.
@@ -102,9 +102,9 @@ Particle Accelerator & Beam Physics
 Scientific works in particle and beam modeling.
 
 #. Sandberg R T, Lehe R, Mitchell C E, Garten M, Myers A, Qiang J, Vay J-L and Huebl A.
-  **Synthesizing Particle-in-Cell Simulations Through Learning and GPU Computing for Hybrid Particle Accelerator Beamlines**.
-  Proc. of Platform for Advanced Scientific Computing (PASC'24), *submitted*, 2024.
-  `preprint <http://arxiv.org/abs/2402.17248>__`
+   **Synthesizing Particle-in-Cell Simulations Through Learning and GPU Computing for Hybrid Particle Accelerator Beamlines**.
+   Proc. of Platform for Advanced Scientific Computing (PASC'24), *submitted*, 2024.
+   `preprint <http://arxiv.org/abs/2402.17248>`__
 
 #. Sandberg R T, Lehe R, Mitchell C E, Garten M, Qiang J, Vay J-L, Huebl A.
    **Hybrid Beamline Element ML-Training for Surrogates in the ImpactX Beam-Dynamics Code**.
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 5f97730e880..40c16a06df5 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -2689,10 +2689,11 @@ In-situ capabilities can be used by turning on Sensei or Ascent (provided they a
 
 * ``<diag_name>.<species_name>.variables`` (list of `strings` separated by spaces, optional)
     List of particle quantities to write to output.
-    Choices are ``w`` for the particle weight and ``ux`` ``uy`` ``uz`` for the particle momenta.
+    Choices are ``x``, ``y``, ``z`` for the particle positions (3D and RZ), ``x`` & ``z`` in 2D, ``z`` in 1D,
+    ``w`` for the particle weight and ``ux``, ``uy``, ``uz`` for the particle momenta.
     When using the lab-frame electrostatic solver, ``phi`` (electrostatic potential, on the macroparticles) is also available.
     By default, all particle quantities (except ``phi``) are written.
-    If ``<diag_name>.<species_name>.variables = none``, no particle data are written, except for particle positions, which are always included.
+    If ``<diag_name>.<species_name>.variables = none``, no particle data are written.
 
 * ``<diag_name>.<species_name>.random_fraction`` (`float`) optional
     If provided ``<diag_name>.<species_name>.random_fraction = a``, only `a` fraction of the particle data of this species will be dumped randomly in diag ``<diag_name>``, i.e. if `rand() < a`, this particle will be dumped, where `rand()` denotes a random number generator.
diff --git a/Examples/Physics_applications/laser_acceleration/inputs_rz b/Examples/Physics_applications/laser_acceleration/inputs_rz
index 7365523b574..b61ed6f83ec 100644
--- a/Examples/Physics_applications/laser_acceleration/inputs_rz
+++ b/Examples/Physics_applications/laser_acceleration/inputs_rz
@@ -89,5 +89,5 @@ diagnostics.diags_names = diag1
 diag1.intervals = 10
 diag1.diag_type = Full
 diag1.fields_to_plot = Er Et Ez Br Bt Bz jr jt jz rho
-diag1.electrons.variables = w ux uy uz orig_x orig_z
-diag1.beam.variables = w ux uy uz
+diag1.electrons.variables = x y z w ux uy uz orig_x orig_z
+diag1.beam.variables = x y z w ux uy uz
diff --git a/Examples/Physics_applications/uniform_plasma/inputs_3d b/Examples/Physics_applications/uniform_plasma/inputs_3d
index b0275f6eed6..83310088176 100644
--- a/Examples/Physics_applications/uniform_plasma/inputs_3d
+++ b/Examples/Physics_applications/uniform_plasma/inputs_3d
@@ -45,7 +45,7 @@ electrons.uz_th  = 0.01 # uth the std of the (unitless) momentum
 diagnostics.diags_names = diag1 chk
 diag1.intervals = 4
 diag1.diag_type = Full
-diag1.electrons.variables = ux uy uz w
+diag1.electrons.variables = x y z ux uy uz w
 diag1.fields_to_plot = Bx By Bz Ex Ey Ez jx jy jz rho
 
 chk.intervals = 6
diff --git a/Examples/Tests/Implicit/inputs_1d b/Examples/Tests/Implicit/inputs_1d
index 50d28a2db75..465d9dd6965 100644
--- a/Examples/Tests/Implicit/inputs_1d
+++ b/Examples/Tests/Implicit/inputs_1d
@@ -71,8 +71,8 @@ diagnostics.diags_names = diag1
 diag1.intervals = 100
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho divE
-diag1.electrons.variables = w ux uy uz
-diag1.protons.variables = w ux uy uz
+diag1.electrons.variables = z w ux uy uz
+diag1.protons.variables = z w ux uy uz
 
 warpx.reduced_diags_names = particle_energy field_energy
 particle_energy.type = ParticleEnergy
diff --git a/Examples/Tests/Implicit/inputs_1d_semiimplicit b/Examples/Tests/Implicit/inputs_1d_semiimplicit
index 2271a0bb1bc..4008a559588 100644
--- a/Examples/Tests/Implicit/inputs_1d_semiimplicit
+++ b/Examples/Tests/Implicit/inputs_1d_semiimplicit
@@ -71,8 +71,8 @@ diagnostics.diags_names = diag1
 diag1.intervals = 100
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho divE
-diag1.electrons.variables = w ux uy uz
-diag1.protons.variables = w ux uy uz
+diag1.electrons.variables = z w ux uy uz
+diag1.protons.variables = z w ux uy uz
 
 warpx.reduced_diags_names = particle_energy field_energy
 particle_energy.type = ParticleEnergy
diff --git a/Examples/Tests/Implicit/inputs_vandb_2d b/Examples/Tests/Implicit/inputs_vandb_2d
index 2dc57323efe..33ce964710a 100644
--- a/Examples/Tests/Implicit/inputs_vandb_2d
+++ b/Examples/Tests/Implicit/inputs_vandb_2d
@@ -82,8 +82,8 @@ diagnostics.diags_names = diag1
 diag1.intervals = 20
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho divE
-diag1.electrons.variables = w ux uy uz
-diag1.protons.variables = w ux uy uz
+diag1.electrons.variables = x z w ux uy uz
+diag1.protons.variables = x z w ux uy uz
 
 warpx.reduced_diags_names = particle_energy field_energy
 particle_energy.type = ParticleEnergy
diff --git a/Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py b/Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py
index 8e05bfa0c68..e63a3f02ba8 100755
--- a/Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py
+++ b/Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py
@@ -18,6 +18,7 @@
 inverse t(r) can be solved for exactly.
 """
 import os
+import re
 import sys
 
 import numpy as np
@@ -36,6 +37,15 @@
 ds = yt.load( filename )
 t_max = ds.current_time.item()  # time of simulation
 
+# Parse test name and check if particle_shape = 4 is used
+emass_10 = True if re.search('emass_10', filename) else False
+
+if emass_10:
+    l2_tolerance = 0.096
+    m_e = 10
+else:
+    l2_tolerance = 0.05
+    m_e = 9.10938356e-31 #Electron mass in kg
 ndims = np.count_nonzero(ds.domain_dimensions > 1)
 
 if ndims == 2:
@@ -60,7 +70,6 @@
 
 # Constants
 eps_0 = 8.8541878128e-12  #Vacuum Permittivity in C/(V*m)
-m_e = 9.10938356e-31  #Electron mass in kg
 q_e = -1.60217662e-19  #Electron charge in C
 pi = np.pi  #Circular constant of the universe
 r_0 = 0.1  #Initial radius of sphere
@@ -134,9 +143,9 @@ def calculate_error(E_axis, xmin, dx, nx):
 print("L2 error along y-axis = %s" %L2_error_y)
 print("L2 error along z-axis = %s" %L2_error_z)
 
-assert L2_error_x < 0.05
-assert L2_error_y < 0.05
-assert L2_error_z < 0.05
+assert L2_error_x < l2_tolerance
+assert L2_error_y < l2_tolerance
+assert L2_error_z < l2_tolerance
 
 # Check conservation of energy
 def return_energies(iteration):
diff --git a/Examples/Tests/electrostatic_sphere/inputs_rz b/Examples/Tests/electrostatic_sphere/inputs_rz
index 26d8b0c7afc..2b6151e6d8c 100644
--- a/Examples/Tests/electrostatic_sphere/inputs_rz
+++ b/Examples/Tests/electrostatic_sphere/inputs_rz
@@ -39,5 +39,5 @@ diag1.fields_to_plot = Er Et Ez rho
 diag2.intervals = 30
 diag2.diag_type = Full
 diag2.fields_to_plot = none
-diag2.electron.variables = ux uy uz w phi
+diag2.electron.variables = x y z ux uy uz w phi
 diag2.format = openpmd
diff --git a/Examples/Tests/langmuir/inputs_3d b/Examples/Tests/langmuir/inputs_3d
index 432b4788171..604ed9bc514 100644
--- a/Examples/Tests/langmuir/inputs_3d
+++ b/Examples/Tests/langmuir/inputs_3d
@@ -93,5 +93,5 @@ diagnostics.diags_names = diag1
 diag1.intervals = max_step
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell rho
-diag1.electrons.variables = w ux
-diag1.positrons.variables = uz
+diag1.electrons.variables = x y z w ux
+diag1.positrons.variables = x y z uz
diff --git a/Examples/Tests/photon_pusher/inputs_3d b/Examples/Tests/photon_pusher/inputs_3d
index e2c75981e0e..3977a187e1b 100644
--- a/Examples/Tests/photon_pusher/inputs_3d
+++ b/Examples/Tests/photon_pusher/inputs_3d
@@ -140,19 +140,19 @@ p_dn_10.single_particle_weight = 1.0
 diagnostics.diags_names = diag1
 diag1.intervals = 50
 diag1.diag_type = Full
-diag1.p_xp_1.variables = ux uy uz
-diag1.p_xn_1.variables = ux uy uz
-diag1.p_yp_1.variables = ux uy uz
-diag1.p_yn_1.variables = ux uy uz
-diag1.p_zp_1.variables = ux uy uz
-diag1.p_zn_1.variables = ux uy uz
-diag1.p_dp_1.variables = ux uy uz
-diag1.p_dn_1.variables = ux uy uz
-diag1.p_xp_10.variables = ux uy uz
-diag1.p_xn_10.variables = ux uy uz
-diag1.p_yp_10.variables = ux uy uz
-diag1.p_yn_10.variables = ux uy uz
-diag1.p_zp_10.variables = ux uy uz
-diag1.p_zn_10.variables = ux uy uz
-diag1.p_dp_10.variables = ux uy uz
-diag1.p_dn_10.variables = ux uy uz
+diag1.p_xp_1.variables = x y z ux uy uz
+diag1.p_xn_1.variables = x y z ux uy uz
+diag1.p_yp_1.variables = x y z ux uy uz
+diag1.p_yn_1.variables = x y z ux uy uz
+diag1.p_zp_1.variables = x y z ux uy uz
+diag1.p_zn_1.variables = x y z ux uy uz
+diag1.p_dp_1.variables = x y z ux uy uz
+diag1.p_dn_1.variables = x y z ux uy uz
+diag1.p_xp_10.variables = x y z ux uy uz
+diag1.p_xn_10.variables = x y z ux uy uz
+diag1.p_yp_10.variables = x y z ux uy uz
+diag1.p_yn_10.variables = x y z ux uy uz
+diag1.p_zp_10.variables = x y z ux uy uz
+diag1.p_zn_10.variables = x y z ux uy uz
+diag1.p_dp_10.variables = x y z ux uy uz
+diag1.p_dn_10.variables = x y z ux uy uz
diff --git a/Examples/Tests/plasma_lens/inputs_3d b/Examples/Tests/plasma_lens/inputs_3d
index 7763013adb1..a87c1eb7f44 100644
--- a/Examples/Tests/plasma_lens/inputs_3d
+++ b/Examples/Tests/plasma_lens/inputs_3d
@@ -48,4 +48,4 @@ particles.repeated_plasma_lens_strengths_B = 0.0 0.0 0.0 0.0
 diagnostics.diags_names = diag1
 diag1.intervals = 84
 diag1.diag_type = Full
-diag1.electrons.variables = ux uy uz
+diag1.electrons.variables = x y z ux uy uz
diff --git a/Examples/Tests/plasma_lens/inputs_boosted_3d b/Examples/Tests/plasma_lens/inputs_boosted_3d
index ecdf3b3baec..fa18ac439c4 100644
--- a/Examples/Tests/plasma_lens/inputs_boosted_3d
+++ b/Examples/Tests/plasma_lens/inputs_boosted_3d
@@ -50,4 +50,4 @@ particles.repeated_plasma_lens_strengths_B = 0.0 0.0 0.0 0.0
 diagnostics.diags_names = diag1
 diag1.intervals = 84
 diag1.diag_type = Full
-diag1.electrons.variables = ux uy uz
+diag1.electrons.variables = x y z ux uy uz
diff --git a/Examples/Tests/plasma_lens/inputs_lattice_3d b/Examples/Tests/plasma_lens/inputs_lattice_3d
index a41ca8baee9..94e57de8290 100644
--- a/Examples/Tests/plasma_lens/inputs_lattice_3d
+++ b/Examples/Tests/plasma_lens/inputs_lattice_3d
@@ -70,4 +70,4 @@ plasmalens4.dEdx = 200000.
 diagnostics.diags_names = diag1
 diag1.intervals = 84
 diag1.diag_type = Full
-diag1.electrons.variables = ux uy uz
+diag1.electrons.variables = x y z ux uy uz
diff --git a/Examples/Tests/plasma_lens/inputs_short_3d b/Examples/Tests/plasma_lens/inputs_short_3d
index 8fa111df3c9..5c91ab26498 100644
--- a/Examples/Tests/plasma_lens/inputs_short_3d
+++ b/Examples/Tests/plasma_lens/inputs_short_3d
@@ -52,4 +52,4 @@ particles.repeated_plasma_lens_strengths_B = 0.0 0.0 0.0 0.0
 diagnostics.diags_names = diag1
 diag1.intervals = 1 #84
 diag1.diag_type = Full
-diag1.electrons.variables = ux uy uz
+diag1.electrons.variables = x y z ux uy uz
diff --git a/Examples/Tests/qed/breit_wheeler/inputs_2d b/Examples/Tests/qed/breit_wheeler/inputs_2d
index 7f507dd08d7..857b3243ac6 100644
--- a/Examples/Tests/qed/breit_wheeler/inputs_2d
+++ b/Examples/Tests/qed/breit_wheeler/inputs_2d
@@ -201,19 +201,19 @@ diagnostics.diags_names = diag1
 diag1.intervals = 2
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex
-diag1.p1.variables = ux uy uz w opticalDepthBW
-diag1.p2.variables = ux uy uz w opticalDepthBW
-diag1.p3.variables = ux uy uz w opticalDepthBW
-diag1.p4.variables = ux uy uz w opticalDepthBW
-
-diag1.ele1.variables = ux uy uz w opticalDepthQSR
-diag1.ele2.variables = ux uy uz w opticalDepthQSR
-diag1.ele3.variables = ux uy uz w opticalDepthQSR
-diag1.ele4.variables = ux uy uz w opticalDepthQSR
-
-diag1.pos1.variables = ux uy uz w opticalDepthQSR
-diag1.pos2.variables = ux uy uz w opticalDepthQSR
-diag1.pos3.variables = ux uy uz w opticalDepthQSR
-diag1.pos4.variables = ux uy uz w opticalDepthQSR
+diag1.p1.variables = x z ux uy uz w opticalDepthBW
+diag1.p2.variables = x z ux uy uz w opticalDepthBW
+diag1.p3.variables = x z ux uy uz w opticalDepthBW
+diag1.p4.variables = x z ux uy uz w opticalDepthBW
+
+diag1.ele1.variables = x z ux uy uz w opticalDepthQSR
+diag1.ele2.variables = x z ux uy uz w opticalDepthQSR
+diag1.ele3.variables = x z ux uy uz w opticalDepthQSR
+diag1.ele4.variables = x z ux uy uz w opticalDepthQSR
+
+diag1.pos1.variables = x z ux uy uz w opticalDepthQSR
+diag1.pos2.variables = x z ux uy uz w opticalDepthQSR
+diag1.pos3.variables = x z ux uy uz w opticalDepthQSR
+diag1.pos4.variables = x z ux uy uz w opticalDepthQSR
 
 diag1.format = plotfile
diff --git a/Examples/Tests/qed/breit_wheeler/inputs_3d b/Examples/Tests/qed/breit_wheeler/inputs_3d
index 2749a60b888..b79dfd9c030 100644
--- a/Examples/Tests/qed/breit_wheeler/inputs_3d
+++ b/Examples/Tests/qed/breit_wheeler/inputs_3d
@@ -201,19 +201,19 @@ diagnostics.diags_names = diag1
 diag1.intervals = 2
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex
-diag1.p1.variables = ux uy uz w opticalDepthBW
-diag1.p2.variables = ux uy uz w opticalDepthBW
-diag1.p3.variables = ux uy uz w opticalDepthBW
-diag1.p4.variables = ux uy uz w opticalDepthBW
-
-diag1.ele1.variables = ux uy uz w opticalDepthQSR
-diag1.ele2.variables = ux uy uz w opticalDepthQSR
-diag1.ele3.variables = ux uy uz w opticalDepthQSR
-diag1.ele4.variables = ux uy uz w opticalDepthQSR
-
-diag1.pos1.variables = ux uy uz w opticalDepthQSR
-diag1.pos2.variables = ux uy uz w opticalDepthQSR
-diag1.pos3.variables = ux uy uz w opticalDepthQSR
-diag1.pos4.variables = ux uy uz w opticalDepthQSR
+diag1.p1.variables = x y z ux uy uz w opticalDepthBW
+diag1.p2.variables = x y z ux uy uz w opticalDepthBW
+diag1.p3.variables = x y z ux uy uz w opticalDepthBW
+diag1.p4.variables = x y z ux uy uz w opticalDepthBW
+
+diag1.ele1.variables = x y z ux uy uz w opticalDepthQSR
+diag1.ele2.variables = x y z ux uy uz w opticalDepthQSR
+diag1.ele3.variables = x y z ux uy uz w opticalDepthQSR
+diag1.ele4.variables = x y z ux uy uz w opticalDepthQSR
+
+diag1.pos1.variables = x y z ux uy uz w opticalDepthQSR
+diag1.pos2.variables = x y z ux uy uz w opticalDepthQSR
+diag1.pos3.variables = x y z ux uy uz w opticalDepthQSR
+diag1.pos4.variables = x y z ux uy uz w opticalDepthQSR
 
 diag1.format = plotfile
diff --git a/Examples/Tests/qed/quantum_synchrotron/inputs_2d b/Examples/Tests/qed/quantum_synchrotron/inputs_2d
index 2fc00231280..2ac2c782ccd 100644
--- a/Examples/Tests/qed/quantum_synchrotron/inputs_2d
+++ b/Examples/Tests/qed/quantum_synchrotron/inputs_2d
@@ -178,12 +178,12 @@ diagnostics.diags_names = diag1
 diag1.intervals = 2
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex
-diag1.p1.variables = ux uy uz w opticalDepthQSR
-diag1.p2.variables = ux uy uz w opticalDepthQSR
-diag1.p3.variables = ux uy uz w opticalDepthQSR
-diag1.p4.variables = ux uy uz w opticalDepthQSR
-
-diag1.qsp_1.variables = ux uy uz w opticalDepthBW
-diag1.qsp_2.variables = ux uy uz w opticalDepthBW
-diag1.qsp_3.variables = ux uy uz w opticalDepthBW
-diag1.qsp_4.variables = ux uy uz w opticalDepthBW
+diag1.p1.variables = x z ux uy uz w opticalDepthQSR
+diag1.p2.variables = x z ux uy uz w opticalDepthQSR
+diag1.p3.variables = x z ux uy uz w opticalDepthQSR
+diag1.p4.variables = x z ux uy uz w opticalDepthQSR
+
+diag1.qsp_1.variables = x z ux uy uz w opticalDepthBW
+diag1.qsp_2.variables = x z ux uy uz w opticalDepthBW
+diag1.qsp_3.variables = x z ux uy uz w opticalDepthBW
+diag1.qsp_4.variables = x z ux uy uz w opticalDepthBW
diff --git a/Examples/Tests/qed/quantum_synchrotron/inputs_3d b/Examples/Tests/qed/quantum_synchrotron/inputs_3d
index 9bc963f7e94..429666ef938 100644
--- a/Examples/Tests/qed/quantum_synchrotron/inputs_3d
+++ b/Examples/Tests/qed/quantum_synchrotron/inputs_3d
@@ -178,12 +178,12 @@ diagnostics.diags_names = diag1
 diag1.intervals = 2
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex
-diag1.p1.variables = ux uy uz w opticalDepthQSR
-diag1.p2.variables = ux uy uz w opticalDepthQSR
-diag1.p3.variables = ux uy uz w opticalDepthQSR
-diag1.p4.variables = ux uy uz w opticalDepthQSR
-
-diag1.qsp_1.variables = ux uy uz w opticalDepthBW
-diag1.qsp_2.variables = ux uy uz w opticalDepthBW
-diag1.qsp_3.variables = ux uy uz w opticalDepthBW
-diag1.qsp_4.variables = ux uy uz w opticalDepthBW
+diag1.p1.variables = x y z ux uy uz w opticalDepthQSR
+diag1.p2.variables = x y z ux uy uz w opticalDepthQSR
+diag1.p3.variables = x y z ux uy uz w opticalDepthQSR
+diag1.p4.variables = x y z ux uy uz w opticalDepthQSR
+
+diag1.qsp_1.variables = x y z ux uy uz w opticalDepthBW
+diag1.qsp_2.variables = x y z ux uy uz w opticalDepthBW
+diag1.qsp_3.variables = x y z ux uy uz w opticalDepthBW
+diag1.qsp_4.variables = x y z ux uy uz w opticalDepthBW
diff --git a/Examples/Tests/resampling/inputs_1d_velocity_coincidence_thinning_cartesian b/Examples/Tests/resampling/inputs_1d_velocity_coincidence_thinning_cartesian
new file mode 100644
index 00000000000..97830d8bb9e
--- /dev/null
+++ b/Examples/Tests/resampling/inputs_1d_velocity_coincidence_thinning_cartesian
@@ -0,0 +1,45 @@
+max_step = 4
+warpx.verbose = 1
+warpx.const_dt = 1e-10
+amr.n_cell = 256
+amr.max_grid_size = 64
+amr.max_level = 0
+geometry.dims = 1
+geometry.prob_lo = 0
+geometry.prob_hi = 0.1
+
+# Boundary condition and field solver
+boundary.field_lo = periodic
+boundary.field_hi = periodic
+boundary.particle_lo = periodic
+boundary.particle_hi = periodic
+algo.particle_shape = 1
+algo.maxwell_solver = none
+
+particles.species_names = hydrogen
+hydrogen.mass = 1.67262192369e-27
+hydrogen.charge = 1.602176634e-19
+hydrogen.injection_style = nrandompercell
+hydrogen.initialize_self_fields = 0
+hydrogen.do_not_push = 1
+hydrogen.do_resampling = 1
+hydrogen.resampling_min_ppc = 10
+hydrogen.resampling_trigger_intervals = 1::2
+hydrogen.resampling_algorithm = velocity_coincidence_thinning
+hydrogen.resampling_algorithm_velocity_grid_type = cartesian
+hydrogen.resampling_algorithm_delta_u = 10000000.0 5000000.0 7000000.0
+hydrogen.num_particles_per_cell = 2500
+hydrogen.momentum_distribution_type = gaussian
+hydrogen.ux_m = 0.0
+hydrogen.uy_m = 0.0
+hydrogen.uz_m = 0.0
+hydrogen.ux_th = 0.000326
+hydrogen.uy_th = 0.000326
+hydrogen.uz_th = 0.000326
+hydrogen.profile = constant
+hydrogen.density = 1e+19
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.intervals = 4
+diag1.diag_type = Full
diff --git a/NOTICE.txt b/NOTICE.txt
new file mode 100644
index 00000000000..e556dddec6c
--- /dev/null
+++ b/NOTICE.txt
@@ -0,0 +1,19 @@
+This Software was developed under funding from the U.S. Department
+of Energy and the U.S. Government consequently retains certain rights. As
+such, the U.S. Government has been granted for itself and others acting on
+its behalf a paid-up, nonexclusive, irrevocable, worldwide license in the
+Software to reproduce, distribute copies to the public, prepare derivative
+works, and perform publicly and display publicly, and to permit other to do
+so.
+
+Reference herein to any specific commercial product, process, or service
+by trade name, trademark, manufacturer, or otherwise does not necessarily
+constitute or imply its endorsement, recommendation, or favoring by the
+United States Government, The Regents of the University of California, or
+Lawrence Berkeley National Laboratory.
+
+The views and opinions of authors expressed herein do not necessarily
+state or reflect those of the United States Government,
+The Regents of the University of California, or
+Lawrence Berkeley National Laboratory, and shall not be used for advertising
+or product endorsement purposes.
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index 03f2ed14732..e7575f83b26 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -254,7 +254,7 @@ def init(self, kw):
         self.resampling_algorithm_n_theta = kw.pop('warpx_resampling_algorithm_n_theta', None)
         self.resampling_algorithm_n_phi = kw.pop('warpx_resampling_algorithm_n_phi', None)
         self.resampling_algorithm_delta_u = kw.pop('warpx_resampling_algorithm_delta_u', None)
-        if np.size(self.resampling_algorithm_delta_u) == 1:
+        if self.resampling_algorithm_delta_u is not None and np.size(self.resampling_algorithm_delta_u) == 1:
             self.resampling_algorithm_delta_u = [self.resampling_algorithm_delta_u]*3
 
     def species_initialize_inputs(self, layout,
diff --git a/README.md b/README.md
index 87cf2b012c0..b2de2211b60 100644
--- a/README.md
+++ b/README.md
@@ -56,13 +56,6 @@ If you have questions about your rights to use or distribute this software,
 please contact Berkeley Lab's Innovation & Partnerships Office at
 IPO@lbl.gov.
 
-NOTICE.  This Software was developed under funding from the U.S. Department
-of Energy and the U.S. Government consequently retains certain rights. As
-such, the U.S. Government has been granted for itself and others acting on
-its behalf a paid-up, nonexclusive, irrevocable, worldwide license in the
-Software to reproduce, distribute copies to the public, prepare derivative
-works, and perform publicly and display publicly, and to permit other to do
-so.
-
 Please see the full license agreement in [LICENSE.txt](LICENSE.txt).
+Please see the notices in [NOTICE.txt](NOTICE.txt).
 The SPDX license identifier is `BSD-3-Clause-LBNL`.
diff --git a/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame_MR_emass_10.json b/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame_MR_emass_10.json
new file mode 100644
index 00000000000..024127a1bd2
--- /dev/null
+++ b/Regression/Checksum/benchmarks_json/ElectrostaticSphereLabFrame_MR_emass_10.json
@@ -0,0 +1,23 @@
+{
+  "lev=0": {
+    "Ex": 7.110909624093144,
+    "Ey": 7.110909624093145,
+    "Ez": 7.110909624093143,
+    "rho": 0.0
+  },
+  "lev=1": {
+    "Ex": 14.281015560380963,
+    "Ey": 14.281015560380965,
+    "Ez": 14.281015560380965,
+    "rho": 2.6092568008333786e-10
+  },
+  "electron": {
+    "particle_momentum_x": 1.80842228672388e-24,
+    "particle_momentum_y": 1.8084222867238806e-24,
+    "particle_momentum_z": 1.7598771525647628e-24,
+    "particle_position_x": 327.46875,
+    "particle_position_y": 327.46875,
+    "particle_position_z": 327.46875,
+    "particle_weight": 6212.501525878906
+  }
+}
\ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/resample_velocity_coincidence_thinning_cartesian.json b/Regression/Checksum/benchmarks_json/resample_velocity_coincidence_thinning_cartesian.json
new file mode 100644
index 00000000000..9f1dae16e7f
--- /dev/null
+++ b/Regression/Checksum/benchmarks_json/resample_velocity_coincidence_thinning_cartesian.json
@@ -0,0 +1,20 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.0,
+    "Bz": 0.0,
+    "Ex": 0.0,
+    "Ey": 0.0,
+    "Ez": 0.0,
+    "jx": 434240.88525516074,
+    "jy": 315768.9402331213,
+    "jz": 642846.5589763735
+  },
+  "hydrogen": {
+    "particle_momentum_x": 6.919253235905925e-20,
+    "particle_momentum_y": 7.069867215063566e-20,
+    "particle_momentum_z": 7.606884436016626e-20,
+    "particle_position_x": 25.60006676227001,
+    "particle_weight": 1.0000000000000001e+18
+  }
+}
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index f8859ef7df5..0f3957a8cce 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -465,6 +465,23 @@ doVis = 0
 compareParticles = 0
 analysisRoutine = Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py
 
+[ElectrostaticSphereLabFrame_MR_emass_10]
+buildDir = .
+inputFile = Examples/Tests/electrostatic_sphere/inputs_3d
+runtime_params = warpx.do_electrostatic=labframe diag2.electron.variables=x y z ux uy uz w warpx.abort_on_warning_threshold=medium electron.mass = 10 amr.max_level = 1 amr.ref_ratio_vect = 2 2 2 warpx.fine_tag_lo = -0.5 -0.5 -0.5 warpx.fine_tag_hi =  0.5  0.5  0.5 max_step = 2
+dim = 3
+addToCompileString =
+cmakeSetupOpts = -DWarpX_DIMS=3
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/electrostatic_sphere/analysis_electrostatic_sphere.py
+
 [ElectrostaticSphereEB]
 buildDir = .
 inputFile = Examples/Tests/electrostatic_sphere_eb/inputs_3d
@@ -554,7 +571,7 @@ analysisRoutine = Examples/Tests/electrostatic_sphere_eb/analysis_rz.py
 [ElectrostaticSphereLabFrame]
 buildDir = .
 inputFile = Examples/Tests/electrostatic_sphere/inputs_3d
-runtime_params = warpx.do_electrostatic=labframe diag2.electron.variables=ux uy uz w phi
+runtime_params = warpx.do_electrostatic=labframe diag2.electron.variables=x y z ux uy uz w phi
 dim = 3
 addToCompileString =
 cmakeSetupOpts = -DWarpX_DIMS=3
@@ -1189,7 +1206,7 @@ analysisOutputImage = langmuir_fluid_multi_analysis.png
 [Langmuir_multi_1d]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_1d
-runtime_params = algo.current_deposition=esirkepov diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz
+runtime_params = algo.current_deposition=esirkepov diag1.electrons.variables=z w ux uy uz diag1.positrons.variables=z w ux uy uz
 dim = 1
 addToCompileString = USE_OPENPMD=TRUE QED=FALSE
 cmakeSetupOpts = -DWarpX_DIMS=1 -DWarpX_OPENPMD=ON -DWarpX_QED=OFF
@@ -1208,7 +1225,7 @@ analysisOutputImage = langmuir_multi_1d_analysis.png
 [Langmuir_multi_2d_MR]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver = ckc  warpx.use_filter = 1  amr.max_level = 1  amr.ref_ratio = 4  warpx.fine_tag_lo = -10.e-6 -10.e-6  warpx.fine_tag_hi = 10.e-6 10.e-6  diag1.electrons.variables = w ux uy uz  diag1.positrons.variables = w ux uy uz
+runtime_params = algo.maxwell_solver = ckc  warpx.use_filter = 1  amr.max_level = 1  amr.ref_ratio = 4  warpx.fine_tag_lo = -10.e-6 -10.e-6  warpx.fine_tag_hi = 10.e-6 10.e-6  diag1.electrons.variables = x z w ux uy uz  diag1.positrons.variables = x z w ux uy uz
 dim = 2
 addToCompileString =
 cmakeSetupOpts = -DWarpX_DIMS=2
@@ -1227,7 +1244,7 @@ analysisOutputImage = Langmuir_multi_2d_MR.png
 [Langmuir_multi_2d_MR_anisotropic]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver = ckc  warpx.use_filter = 1  amr.max_level = 1 amr.ref_ratio_vect = 4 2 warpx.fine_tag_lo = -10.e-6 -10.e-6  warpx.fine_tag_hi = 10.e-6 10.e-6  diag1.electrons.variables = w ux uy uz  diag1.positrons.variables = w ux uy uz
+runtime_params = algo.maxwell_solver = ckc  warpx.use_filter = 1  amr.max_level = 1 amr.ref_ratio_vect = 4 2 warpx.fine_tag_lo = -10.e-6 -10.e-6  warpx.fine_tag_hi = 10.e-6 10.e-6  diag1.electrons.variables = x z w ux uy uz  diag1.positrons.variables = x z w ux uy uz
 dim = 2
 addToCompileString =
 cmakeSetupOpts = -DWarpX_DIMS=2
@@ -1246,7 +1263,7 @@ analysisOutputImage = Langmuir_multi_2d_MR.png
 [Langmuir_multi_2d_MR_momentum_conserving]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver=ckc warpx.use_filter=1 amr.max_level=1 amr.ref_ratio=4 warpx.fine_tag_lo=-10.e-6 -10.e-6 warpx.fine_tag_hi=10.e-6 10.e-6 algo.field_gathering=momentum-conserving diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz
+runtime_params = algo.maxwell_solver=ckc warpx.use_filter=1 amr.max_level=1 amr.ref_ratio=4 warpx.fine_tag_lo=-10.e-6 -10.e-6 warpx.fine_tag_hi=10.e-6 10.e-6 algo.field_gathering=momentum-conserving diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz
 dim = 2
 addToCompileString =
 cmakeSetupOpts = -DWarpX_DIMS=2
@@ -1265,7 +1282,7 @@ analysisOutputImage = Langmuir_multi_2d_MR_momentum_conserving.png
 [Langmuir_multi_2d_MR_psatd]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver = psatd  warpx.use_filter = 1  amr.max_level = 1  amr.ref_ratio = 4  warpx.fine_tag_lo = -10.e-6 -10.e-6  warpx.fine_tag_hi = 10.e-6 10.e-6  diag1.electrons.variables = w ux uy uz  diag1.positrons.variables = w ux uy uz psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
+runtime_params = algo.maxwell_solver = psatd  warpx.use_filter = 1  amr.max_level = 1  amr.ref_ratio = 4  warpx.fine_tag_lo = -10.e-6 -10.e-6  warpx.fine_tag_hi = 10.e-6 10.e-6  diag1.electrons.variables = x z w ux uy uz  diag1.positrons.variables = x z w ux uy uz psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
@@ -1284,7 +1301,7 @@ analysisOutputImage = Langmuir_multi_2d_MR_psatd.png
 [Langmuir_multi_2d_nodal]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = warpx.grid_type=collocated algo.current_deposition=direct diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz
+runtime_params = warpx.grid_type=collocated algo.current_deposition=direct diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz
 dim = 2
 addToCompileString =
 cmakeSetupOpts = -DWarpX_DIMS=2
@@ -1303,7 +1320,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver=psatd diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
+runtime_params = algo.maxwell_solver=psatd diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
@@ -1322,7 +1339,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_current_correction]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=esirkepov psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot =Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=esirkepov psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz diag1.fields_to_plot =Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
@@ -1341,7 +1358,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_current_correction_nodal]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=direct psatd.periodic_single_box_fft=1 psatd.current_correction=1 warpx.grid_type=collocated diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot =Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=direct psatd.periodic_single_box_fft=1 psatd.current_correction=1 warpx.grid_type=collocated diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz diag1.fields_to_plot =Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
@@ -1360,7 +1377,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_momentum_conserving]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver=psatd algo.field_gathering=momentum-conserving diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
+runtime_params = algo.maxwell_solver=psatd algo.field_gathering=momentum-conserving diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
@@ -1417,7 +1434,7 @@ analysisOutputImage = Langmuir_multi_2d_psatd_multiJ_nodal.png
 [Langmuir_multi_2d_psatd_nodal]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver=psatd warpx.grid_type=collocated algo.current_deposition=direct diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
+runtime_params = algo.maxwell_solver=psatd warpx.grid_type=collocated algo.current_deposition=direct diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell warpx.cfl = 0.7071067811865475 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
@@ -1436,7 +1453,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_Vay_deposition]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=vay diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=vay diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
@@ -1455,7 +1472,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_Vay_deposition_particle_shape_4]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=vay diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475 algo.particle_shape=4
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=vay diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475 algo.particle_shape=4
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
@@ -1474,7 +1491,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_Vay_deposition_nodal]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d
-runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 warpx.grid_type=collocated algo.current_deposition=vay diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 warpx.grid_type=collocated algo.current_deposition=vay diag1.electrons.variables=x z w ux uy uz diag1.positrons.variables=x z w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
@@ -1492,7 +1509,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 
 [Langmuir_multi_nodal]
 buildDir = .
-inputFile = Examples/Tests/langmuir/inputs_3d
+inputFile =
 runtime_params = warpx.grid_type=collocated algo.current_deposition=direct
 dim = 3
 addToCompileString =
@@ -1721,7 +1738,7 @@ analysisOutputImage = langmuir_multi_analysis.png
 [Langmuir_multi_rz]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_rz
-runtime_params = diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz diag1.dump_rz_modes=0
+runtime_params = diag1.electrons.variables=x y z w ux uy uz diag1.ions.variables=x y z w ux uy uz diag1.dump_rz_modes=0
 dim = 2
 addToCompileString = USE_RZ=TRUE
 cmakeSetupOpts = -DWarpX_DIMS=RZ
@@ -1741,7 +1758,7 @@ aux1File = Regression/PostProcessingUtils/post_processing_utils.py
 [Langmuir_multi_rz_psatd]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_rz
-runtime_params = algo.maxwell_solver=psatd diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0 psatd.update_with_rho=1 electrons.random_theta=0 ions.random_theta=0 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
+runtime_params = algo.maxwell_solver=psatd diag1.electrons.variables=x y z w ux uy uz diag1.ions.variables=x y z w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0 psatd.update_with_rho=1 electrons.random_theta=0 ions.random_theta=0 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium
 dim = 2
 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE BLAS_LIB=-lblas LAPACK_LIB=-llapack
 cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_PSATD=ON
@@ -1761,7 +1778,7 @@ aux1File = Regression/PostProcessingUtils/post_processing_utils.py
 [Langmuir_multi_rz_psatd_current_correction]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_rz
-runtime_params = algo.maxwell_solver=psatd diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0 amr.max_grid_size=128 psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.fields_to_plot=jr jz Er Ez Bt rho divE electrons.random_theta=0 ions.random_theta=0
+runtime_params = algo.maxwell_solver=psatd diag1.electrons.variables=x y z w ux uy uz diag1.ions.variables=x y z w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0 amr.max_grid_size=128 psatd.periodic_single_box_fft=1 psatd.current_correction=1 diag1.fields_to_plot=jr jz Er Ez Bt rho divE electrons.random_theta=0 ions.random_theta=0
 dim = 2
 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE BLAS_LIB=-lblas LAPACK_LIB=-llapack
 cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_PSATD=ON
@@ -1781,7 +1798,7 @@ aux1File = Regression/PostProcessingUtils/post_processing_utils.py
 [Langmuir_multi_rz_psatd_multiJ]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_rz
-runtime_params = amr.max_grid_size=32 algo.maxwell_solver=psatd diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0 psatd.update_with_rho=1 warpx.n_rz_azimuthal_modes=2 electrons.random_theta=0 electrons.num_particles_per_cell_each_dim=2 4 2 ions.random_theta=0 ions.num_particles_per_cell_each_dim=2 4 2 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=4 warpx.use_filter=1
+runtime_params = amr.max_grid_size=32 algo.maxwell_solver=psatd diag1.electrons.variables=x y z w ux uy uz diag1.ions.variables=x y z w ux uy uz diag1.dump_rz_modes=0 algo.current_deposition=direct warpx.do_dive_cleaning=0 psatd.update_with_rho=1 warpx.n_rz_azimuthal_modes=2 electrons.random_theta=0 electrons.num_particles_per_cell_each_dim=2 4 2 ions.random_theta=0 ions.num_particles_per_cell_each_dim=2 4 2 psatd.current_correction=0 warpx.abort_on_warning_threshold=medium warpx.do_multi_J=1 warpx.do_multi_J_n_depositions=4 warpx.use_filter=1
 dim = 2
 addToCompileString = USE_RZ=TRUE USE_PSATD=TRUE BLAS_LIB=-lblas LAPACK_LIB=-llapack
 cmakeSetupOpts = -DWarpX_DIMS=RZ -DWarpX_PSATD=ON
@@ -4206,6 +4223,23 @@ doVis = 0
 compareParticles = 0
 analysisRoutine = Examples/analysis_default_regression.py
 
+[resample_velocity_coincidence_thinning_cartesian]
+buildDir = .
+inputFile = Examples/Tests/resampling/inputs_1d_velocity_coincidence_thinning_cartesian
+runtime_params =
+dim = 1
+addToCompileString =
+cmakeSetupOpts = -DWarpX_DIMS=1
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/analysis_default_regression.py
+
 [restart]
 buildDir = .
 inputFile = Examples/Tests/restart/inputs
diff --git a/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp b/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp
index 7ed64ccf8ec..1a64ae20f0e 100644
--- a/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp
+++ b/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp
@@ -6,6 +6,8 @@
 #include "Utils/TextMsg.H"
 #include "WarpX.H"
 
+#include <ablastr/warn_manager/WarnManager.H>
+
 #include <AMReX_ParmParse.H>
 
 #include <map>
@@ -13,7 +15,7 @@
 
 using namespace amrex;
 
-ParticleDiag::ParticleDiag(
+ParticleDiag::ParticleDiag (
     const std::string& diag_name, const std::string& name,
     WarpXParticleContainer* pc, PinnedMemoryParticleContainer* pinned_pc):
     m_diag_name(diag_name), m_name(name), m_pc(pc), m_pinned_pc(pinned_pc)
@@ -28,12 +30,17 @@ ParticleDiag::ParticleDiag(
     amrex::Vector<std::string> variables;
     const int variables_specified = pp_diag_name_species_name.queryarr("variables", variables);
 
-    if (variables_specified){
+    if (variables_specified) {
         // If only specific variables have been specified, fill m_plot_flags with zero and only set
         // requested variables to one
         std::fill(m_plot_flags.begin(), m_plot_flags.end(), 0);
+        bool contains_positions = false;
         if (variables[0] != "none"){
-            const std::map<std::string, int> existing_variable_names = pc->getParticleComps();
+            std::map<std::string, int> existing_variable_names = pc->getParticleComps();
+#ifdef WARPX_DIM_RZ
+            // we reconstruct to Cartesian x,y,z for RZ particle output
+            existing_variable_names["y"] = PIdx::theta;
+#endif
             for (const auto& var : variables){
                 if (var == "phi") {
                     // User requests phi on particle. This is *not* part of the variables that
@@ -47,9 +54,21 @@ ParticleDiag::ParticleDiag(
                         "variables argument '" + var
                         +"' is not an existing attribute for this species");
                     m_plot_flags[existing_variable_names.at(var)] = 1;
+
+                    if (var == "x" || var == "y" || var == "z") {
+                        contains_positions = true;
+                    }
                 }
             }
         }
+
+        if (!contains_positions) {
+            ablastr::warn_manager::WMRecordWarning(
+                "Diagnostics",
+                diag_name + "." + name + ".variables contains no particle positions!",
+                ablastr::warn_manager::WarnPriority::high
+            );
+        }
     }
 
 #ifdef WARPX_DIM_RZ
diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp
index c4164e383b1..f1f4d426a4a 100644
--- a/Source/Diagnostics/ParticleIO.cpp
+++ b/Source/Diagnostics/ParticleIO.cpp
@@ -283,7 +283,8 @@ storePhiOnParticles ( PinnedMemoryParticleContainer& tmp,
                     getPosition(ip, xp, yp, zp);
                     int i, j, k;
                     amrex::Real W[AMREX_SPACEDIM][2];
-                    ablastr::particles::compute_weights(xp, yp, zp, plo, dxi, i, j, k, W);
+                    ablastr::particles::compute_weights<amrex::IndexType::NODE>(
+                        xp, yp, zp, plo, dxi, i, j, k, W);
                     amrex::Real const phi_value  = ablastr::particles::interp_field_nodal(i, j, k, W, phi_grid);
                     phi_particle_arr[ip] = phi_value;
                 }
diff --git a/Source/Diagnostics/WarpXOpenPMD.H b/Source/Diagnostics/WarpXOpenPMD.H
index 12aa46de058..99d6e0682ab 100644
--- a/Source/Diagnostics/WarpXOpenPMD.H
+++ b/Source/Diagnostics/WarpXOpenPMD.H
@@ -218,11 +218,13 @@ private:
   /** This function sets up the entries for storing the particle positions and global IDs
   *
   * @param[in] currSpecies Corresponding openPMD species
+  * @param[in] positionComponents user-selected components of the particle position
   * @param[in] np          Number of particles
   * @param[in] isBTD       Is this a back-transformed diagnostics output?
   */
   void SetupPos (
         openPMD::ParticleSpecies& currSpecies,
+        std::vector<std::string> const & positionComponents,
         const unsigned long long& np,
         bool isBTD = false);
 
@@ -231,12 +233,14 @@ private:
    * Sets the entries for storing particle position offset, constant records (charge, mass) and ED-PIC attributes.
    *
    * @param[in] currSpecies Corresponding openPMD species
+   * @param[in] positionComponents user-selected components of the particle position
    * @param[in] np          Number of particles
    * @param[in] charge      Charge of the particles (note: fix for ions)
    * @param[in] mass        Mass of the particles
    */
   void SetConstParticleRecordsEDPIC (
         openPMD::ParticleSpecies& currSpecies,
+        std::vector<std::string> const & positionComponents,
         const unsigned long long& np,
         amrex::ParticleReal charge,
         amrex::ParticleReal mass);
@@ -284,8 +288,8 @@ private:
    * @param[in] name species name
    * @param[in] iteration timestep
    * @param[in] write_real_comp The real attribute ids, from WarpX
-   * @param[in] real_comp_names The real attribute names, from WarpX
    * @param[in] write_int_comp The int attribute ids, from WarpX
+   * @param[in] real_comp_names The real attribute names, from WarpX
    * @param[in] int_comp_names The int attribute names, from WarpX
    * @param[in] charge         Charge of the particles (note: fix for ions)
    * @param[in] mass           Mass of the particles
diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp
index 965eb82b31e..45d5dc1cdee 100644
--- a/Source/Diagnostics/WarpXOpenPMD.cpp
+++ b/Source/Diagnostics/WarpXOpenPMD.cpp
@@ -213,32 +213,29 @@ namespace detail
         return make_pair(record_name, component_name);
     }
 
-    /** Return the component labels for particle positions
+    /** Return the user-selected components for particle positions
+     *
+     * @param[in] write_real_comp The real attribute ids, from WarpX
+     * @param[in] real_comp_names The real attribute names, from WarpX
      */
     inline std::vector< std::string >
-    getParticlePositionComponentLabels (bool ignore_dims=false)
+    getParticlePositionComponentLabels (
+        amrex::Vector<int> const & write_real_comp,
+        amrex::Vector<std::string> const & real_comp_names)
     {
-        using vs = std::vector< std::string >;
-        auto positionComponents = vs{"x", "y", "z"};
-        if (!ignore_dims) {
-#if defined(WARPX_DIM_1D_Z)
-            positionComponents = vs{"z"};
-#elif defined(WARPX_DIM_XZ)
-            positionComponents = vs{"x", "z"};
-#elif defined(WARPX_DIM_RZ)
-            // note: although we internally store particle positions
-            //       for AMReX in r,z and a theta attribute, we
-            //       actually need them for algorithms (e.g. push)
-            //       and I/O in Cartesian.
-            //       Other attributes like momentum are consequently
-            //       stored in x,y,z internally.
-            positionComponents = vs{"x", "y", "z"};
-#elif defined(WARPX_DIM_3D)
-            positionComponents = vs{"x", "y", "z"};
-#else
-#   error Unknown WarpX dimensionality.
-#endif
+        std::vector< std::string > positionComponents;
+
+        int idx = 0;
+        for (auto const & comp : real_comp_names ) {
+            if (write_real_comp[idx]) {
+                if (comp == "position_x" || comp == "position_y" || comp == "position_z") {
+                    std::string const last_letter{comp.back()};
+                    positionComponents.push_back(last_letter);
+                }
+            }
+            idx++;
         }
+
         return positionComponents;
     }
 
@@ -705,16 +702,18 @@ WarpXOpenPMDPlot::DumpToFile (ParticleContainer* pc,
         doParticleSetup = is_first_flush_with_particles || is_last_flush_and_never_particles;
     }
 
+    auto const positionComponents = detail::getParticlePositionComponentLabels(write_real_comp, real_comp_names);
+
     // this setup stage also implicitly calls "makeEmpty" if needed (i.e., is_last_flush_and_never_particles)
     //   for BTD, we call this multiple times as we may resize in subsequent dumps if number of particles in the buffer > 0
     if (doParticleSetup || is_resizing_flush) {
-        SetupPos(currSpecies, NewParticleVectorSize, isBTD);
+        SetupPos(currSpecies, positionComponents, NewParticleVectorSize, isBTD);
         SetupRealProperties(pc, currSpecies, write_real_comp, real_comp_names, write_int_comp, int_comp_names,
                             NewParticleVectorSize, isBTD);
     }
 
     if (is_last_flush_to_step) {
-        SetConstParticleRecordsEDPIC(currSpecies, NewParticleVectorSize, charge, mass);
+        SetConstParticleRecordsEDPIC(currSpecies, positionComponents, NewParticleVectorSize, charge, mass);
     }
 
     // open files from all processors, in case some will not contribute below
@@ -978,6 +977,7 @@ WarpXOpenPMDPlot::SaveRealProperty (ParticleIter& pti,
 void
 WarpXOpenPMDPlot::SetupPos (
     openPMD::ParticleSpecies& currSpecies,
+    std::vector<std::string> const & positionComponents,
     const unsigned long long& np,
     bool const isBTD)
 {
@@ -986,7 +986,6 @@ WarpXOpenPMDPlot::SetupPos (
     auto realType = openPMD::Dataset(openPMD::determineDatatype<amrex::ParticleReal>(), {np}, options);
     auto idType = openPMD::Dataset(openPMD::determineDatatype< uint64_t >(), {np}, options);
 
-    auto const positionComponents = detail::getParticlePositionComponentLabels();
     for( auto const& comp : positionComponents ) {
         currSpecies["position"][comp].resetDataset( realType );
     }
@@ -998,6 +997,7 @@ WarpXOpenPMDPlot::SetupPos (
 void
 WarpXOpenPMDPlot::SetConstParticleRecordsEDPIC (
         openPMD::ParticleSpecies& currSpecies,
+        std::vector<std::string> const & positionComponents,
         const unsigned long long& np,
         amrex::ParticleReal const charge,
         amrex::ParticleReal const mass)
@@ -1006,7 +1006,6 @@ WarpXOpenPMDPlot::SetConstParticleRecordsEDPIC (
     const auto *const scalar = openPMD::RecordComponent::SCALAR;
 
     // define record shape to be number of particles
-    auto const positionComponents = detail::getParticlePositionComponentLabels(true);
     for( auto const& comp : positionComponents ) {
         currSpecies["positionOffset"][comp].resetDataset( realType );
     }
@@ -1015,9 +1014,12 @@ WarpXOpenPMDPlot::SetConstParticleRecordsEDPIC (
 #if defined(WARPX_DIM_1D_Z)
     currSpecies["position"]["x"].resetDataset( realType );
     currSpecies["position"]["y"].resetDataset( realType );
+    currSpecies["positionOffset"]["x"].resetDataset( realType );
+    currSpecies["positionOffset"]["y"].resetDataset( realType );
 #endif
 #if defined(WARPX_DIM_XZ)
     currSpecies["position"]["y"].resetDataset( realType );
+    currSpecies["positionOffset"]["y"].resetDataset( realType );
 #endif
 
     // make constant
@@ -1031,22 +1033,29 @@ WarpXOpenPMDPlot::SetConstParticleRecordsEDPIC (
 #if defined(WARPX_DIM_1D_Z)
     currSpecies["position"]["x"].makeConstant( 0._prt );
     currSpecies["position"]["y"].makeConstant( 0._prt );
+    currSpecies["positionOffset"]["x"].makeConstant( 0._prt );
+    currSpecies["positionOffset"]["y"].makeConstant( 0._prt );
 #endif
 #if defined(WARPX_DIM_XZ)
     currSpecies["position"]["y"].makeConstant( 0._prt );
+    currSpecies["positionOffset"]["y"].makeConstant( 0._prt );
 #endif
 
     // meta data
-    currSpecies["position"].setUnitDimension( detail::getUnitDimension("position") );
-    currSpecies["positionOffset"].setUnitDimension( detail::getUnitDimension("positionOffset") );
+    if (!positionComponents.empty()) {
+        currSpecies["position"].setUnitDimension( detail::getUnitDimension("position") );
+        currSpecies["positionOffset"].setUnitDimension( detail::getUnitDimension("positionOffset") );
+    }
     currSpecies["charge"].setUnitDimension( detail::getUnitDimension("charge") );
     currSpecies["mass"].setUnitDimension( detail::getUnitDimension("mass") );
 
     // meta data for ED-PIC extension
-    currSpecies["position"].setAttribute( "macroWeighted", 0u );
-    currSpecies["position"].setAttribute( "weightingPower", 0.0 );
-    currSpecies["positionOffset"].setAttribute( "macroWeighted", 0u );
-    currSpecies["positionOffset"].setAttribute( "weightingPower", 0.0 );
+    if (!positionComponents.empty()) {
+        currSpecies["position"].setAttribute( "macroWeighted", 0u );
+        currSpecies["position"].setAttribute( "weightingPower", 0.0 );
+        currSpecies["positionOffset"].setAttribute( "macroWeighted", 0u );
+        currSpecies["positionOffset"].setAttribute( "weightingPower", 0.0 );
+    }
     currSpecies["id"].setAttribute( "macroWeighted", 0u );
     currSpecies["id"].setAttribute( "weightingPower", 0.0 );
     currSpecies["charge"].setAttribute( "macroWeighted", 0u );
diff --git a/Source/EmbeddedBoundary/ParticleScraper.H b/Source/EmbeddedBoundary/ParticleScraper.H
index ea4d0bc5c0e..1e915b39381 100644
--- a/Source/EmbeddedBoundary/ParticleScraper.H
+++ b/Source/EmbeddedBoundary/ParticleScraper.H
@@ -183,15 +183,16 @@ scrapeParticlesAtEB (PC& pc, const amrex::Vector<const amrex::MultiFab*>& distan
 
                 int i, j, k;
                 amrex::Real W[AMREX_SPACEDIM][2];
-                ablastr::particles::compute_weights(xp, yp, zp, plo, dxi, i, j, k, W);
+                ablastr::particles::compute_weights<amrex::IndexType::NODE>(
+                    xp, yp, zp, plo, dxi, i, j, k, W);
                 amrex::Real phi_value  = ablastr::particles::interp_field_nodal(i, j, k, W, phi);
 
                 if (phi_value < 0.0)
                 {
                     int ic, jc, kc; // Cell-centered indices
-                    [[maybe_unused]] int nodal;
                     amrex::Real Wc[AMREX_SPACEDIM][2]; // Cell-centered weights
-                    ablastr::particles::compute_weights(xp, yp, zp, plo, dxi, ic, jc, kc, Wc, nodal=0);
+                    ablastr::particles::compute_weights<amrex::IndexType::CELL>(
+                        xp, yp, zp, plo, dxi, ic, jc, kc, Wc);
                     amrex::RealVect normal = DistanceToEB::interp_normal(i, j, k, W, ic, jc, kc, Wc, phi, dxi);
                     DistanceToEB::normalize(normal);
                     amrex::RealVect pos;
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index 51fcb5c262c..de586895c2d 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -645,7 +645,7 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time)
     //    (dt[0] denotes the time step on mesh refinement level 0)
     if (J_in_time == JInTime::Linear)
     {
-        auto& current = (WarpX::do_current_centering) ? current_fp_nodal : current_fp;
+        auto& current = (do_current_centering) ? current_fp_nodal : current_fp;
         mypc->DepositCurrent(current, dt[0], -dt[0]);
         // Synchronize J: filter, exchange boundary, and interpolate across levels.
         // With current centering, the nodal current is deposited in 'current',
@@ -679,7 +679,7 @@ WarpX::OneStep_multiJ (const amrex::Real cur_time)
 
         // Deposit new J at relative time t_deposit_current with time step dt
         // (dt[0] denotes the time step on mesh refinement level 0)
-        auto& current = (WarpX::do_current_centering) ? current_fp_nodal : current_fp;
+        auto& current = (do_current_centering) ? current_fp_nodal : current_fp;
         mypc->DepositCurrent(current, dt[0], t_deposit_current);
         // Synchronize J: filter, exchange boundary, and interpolate across levels.
         // With current centering, the nodal current is deposited in 'current',
diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp
index 5bb50ada278..a74beae3c91 100644
--- a/Source/FieldSolver/ElectrostaticSolver.cpp
+++ b/Source/FieldSolver/ElectrostaticSolver.cpp
@@ -195,7 +195,13 @@ WarpX::AddSpaceChargeField (WarpXParticleContainer& pc)
         pc.DepositCharge(rho, local, reset, apply_boundary_and_scale_volume,
                               interpolate_across_levels);
     }
-
+    for (int lev = 0; lev <= max_level; lev++) {
+        if (lev > 0) {
+            if (charge_buf[lev]) {
+                charge_buf[lev]->setVal(0.);
+            }
+        }
+    }
     SyncRho(rho, rho_coarse, charge_buf); // Apply filter, perform MPI exchange, interpolate across levels
 
     // Get the particle beta vector
@@ -239,7 +245,13 @@ WarpX::AddSpaceChargeFieldLabFrame ()
         int const lev = 0;
         myfl->DepositCharge( lev, *rho_fp[lev] );
     }
-
+    for (int lev = 0; lev <= max_level; lev++) {
+        if (lev > 0) {
+            if (charge_buf[lev]) {
+                charge_buf[lev]->setVal(0.);
+            }
+        }
+    }
     SyncRho(rho_fp, rho_cp, charge_buf); // Apply filter, perform MPI exchange, interpolate across levels
 #ifndef WARPX_DIM_RZ
     for (int lev = 0; lev <= finestLevel(); lev++) {
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index 9daf37af3f8..f3bff69546b 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -882,7 +882,7 @@ WarpX::SyncCurrent (
     WARPX_PROFILE("WarpX::SyncCurrent()");
 
     // If warpx.do_current_centering = 1, center currents from nodal grid to staggered grid
-    if (WarpX::do_current_centering)
+    if (do_current_centering)
     {
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level <= 1,
                                          "warpx.do_current_centering=1 not supported with more than one fine levels");
@@ -1178,7 +1178,7 @@ void WarpX::SumBoundaryJ (
     const amrex::IntVect ng = J.nGrowVect();
     amrex::IntVect ng_depos_J = get_ng_depos_J();
 
-    if (WarpX::do_current_centering)
+    if (do_current_centering)
     {
 #if   defined(WARPX_DIM_1D_Z)
         ng_depos_J[0] += WarpX::current_centering_noz / 2;
diff --git a/Source/Particles/ParticleBoundaryBuffer.cpp b/Source/Particles/ParticleBoundaryBuffer.cpp
index c8c683f0abf..345a0bfa592 100644
--- a/Source/Particles/ParticleBoundaryBuffer.cpp
+++ b/Source/Particles/ParticleBoundaryBuffer.cpp
@@ -98,7 +98,8 @@ struct FindEmbeddedBoundaryIntersection {
                 amrex::Real W[AMREX_SPACEDIM][2];
                 amrex::ParticleReal x_temp=xp, y_temp=yp, z_temp=zp;
                 UpdatePosition(x_temp, y_temp, z_temp, ux, uy, uz, -dt_frac*dt);
-                ablastr::particles::compute_weights(x_temp, y_temp, z_temp, plo, dxi, i, j, k, W);
+                ablastr::particles::compute_weights<amrex::IndexType::NODE>(
+                    x_temp, y_temp, z_temp, plo, dxi, i, j, k, W);
                 amrex::Real phi_value  = ablastr::particles::interp_field_nodal(i, j, k, W, phiarr);
                 return phi_value;
             } );
@@ -115,11 +116,12 @@ struct FindEmbeddedBoundaryIntersection {
         // record the components of the normal on the destination
         int i, j, k;
         amrex::Real W[AMREX_SPACEDIM][2];
-        ablastr::particles::compute_weights(x_temp, y_temp, z_temp, plo, dxi, i, j, k, W);
+        ablastr::particles::compute_weights<amrex::IndexType::NODE>(
+            x_temp, y_temp, z_temp, plo, dxi, i, j, k, W);
         int ic, jc, kc; // Cell-centered indices
-        int nodal;
         amrex::Real Wc[AMREX_SPACEDIM][2]; // Cell-centered weight
-        ablastr::particles::compute_weights(x_temp, y_temp, z_temp, plo, dxi, ic, jc, kc, Wc, nodal=0); // nodal=0 to calculate the weights with respect to the cell-centered nodes
+        ablastr::particles::compute_weights<amrex::IndexType::CELL>(
+            x_temp, y_temp, z_temp, plo, dxi, ic, jc, kc, Wc);
         amrex::RealVect normal = DistanceToEB::interp_normal(i, j, k, W, ic, jc, kc, Wc, phiarr, dxi);
         DistanceToEB::normalize(normal);
 
diff --git a/Source/WarpX.H b/Source/WarpX.H
index c6039387a6e..85afe863662 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -212,7 +212,7 @@ public:
     //! If true, the current is deposited on a nodal grid and then centered onto a staggered grid
     //! using finite centering of order given by #current_centering_nox, #current_centering_noy,
     //! and #current_centering_noz
-    static bool do_current_centering;
+    bool do_current_centering = false;
 
     //! If true, a correction is applied to the current in Fourier space,
     //  to satisfy the continuity equation and charge conservation
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 118654b2b40..40df688a4b4 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -144,8 +144,6 @@ amrex::Vector<FieldBoundaryType> WarpX::field_boundary_hi(AMREX_SPACEDIM,FieldBo
 amrex::Vector<ParticleBoundaryType> WarpX::particle_boundary_lo(AMREX_SPACEDIM,ParticleBoundaryType::Absorbing);
 amrex::Vector<ParticleBoundaryType> WarpX::particle_boundary_hi(AMREX_SPACEDIM,ParticleBoundaryType::Absorbing);
 
-bool WarpX::do_current_centering = false;
-
 int WarpX::n_rz_azimuthal_modes = 1;
 int WarpX::ncomps = 1;
 
diff --git a/Source/ablastr/particles/NodalFieldGather.H b/Source/ablastr/particles/NodalFieldGather.H
index 5a10c73ae20..f324264db7f 100644
--- a/Source/ablastr/particles/NodalFieldGather.H
+++ b/Source/ablastr/particles/NodalFieldGather.H
@@ -12,16 +12,19 @@
 #include <AMReX_Array.H>
 #include <AMReX_Extension.H>
 #include <AMReX_GpuQualifiers.H>
+#include <AMReX_IndexType.H>
 #include <AMReX_Math.H>
 #include <AMReX_REAL.H>
 
+
 namespace ablastr::particles
 {
+
 /**
  * \brief Compute weight of each surrounding node (or cell-centered nodes) in interpolating a nodal ((or a cell-centered node) field
- *        to the given coordinates. If nodal=1, then the calculations will be done with respect to the nodes (default). If nodal=0, then the calculations will be done with respect to the cell-centered nodal)
- *
- * This currently only does linear order.
+ *        to the given coordinates. If template parameter IdxType is amrex::IndexType::CellIndex::NODE, then the calculations will be done with respect to the nodes.
+ *        If template parameter IdxType is amrex::IndexType::CellIndex::CELL, then the calculations will be done with respect to the cell-centered nodal.
+ *        This currently only does linear order.
  *
  * \param xp,yp,zp Particle position coordinates
  * \param plo      Index lower bounds of domain.
@@ -30,24 +33,30 @@ namespace ablastr::particles
  * \param W        2D array of weights to store each neighbouring node (or cell-centered node)
  * \param nodal    Int that tells if the weights are calculated in respect to the nodes (nodal=1) of the cell-centered nodes (nodal=0)
  */
+template <amrex::IndexType::CellIndex IdxType>
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 void compute_weights (const amrex::ParticleReal xp,
                             const amrex::ParticleReal yp,
                             const amrex::ParticleReal zp,
                             amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& plo,
                             amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dxi,
-                            int& i, int& j, int& k, amrex::Real W[AMREX_SPACEDIM][2], int nodal=1) noexcept
+                            int& i, int& j, int& k, amrex::Real W[AMREX_SPACEDIM][2]) noexcept
 {
     using namespace amrex::literals;
 
-#if !((nodal==0)||(nodal==1))
-    ABLASTR_ABORT_WITH_MESSAGE("Error: 'nodal' has to be equal to 0 or 1");
-#endif
+    constexpr auto half_if_cell_centered = [](){
+        if constexpr (IdxType == amrex::IndexType::CellIndex::NODE){
+            return 0.0_rt;
+        }
+        else{
+            return 0.5_rt;
+        }
+    }();
 
 #if (defined WARPX_DIM_3D)
-    const amrex::Real x = (xp - plo[0]) * dxi[0] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
-    const amrex::Real y = (yp - plo[1]) * dxi[1] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
-    const amrex::Real z = (zp - plo[2]) * dxi[2] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
+    const amrex::Real x = (xp - plo[0]) * dxi[0] - half_if_cell_centered;
+    const amrex::Real y = (yp - plo[1]) * dxi[1] - half_if_cell_centered;
+    const amrex::Real z = (zp - plo[2]) * dxi[2] - half_if_cell_centered;
 
     i = static_cast<int>(amrex::Math::floor(x));
     j = static_cast<int>(amrex::Math::floor(y));
@@ -64,17 +73,17 @@ void compute_weights (const amrex::ParticleReal xp,
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
 
 #   if (defined WARPX_DIM_XZ)
-    const amrex::Real x = (xp - plo[0]) * dxi[0] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
+    const amrex::Real x = (xp - plo[0]) * dxi[0] - half_if_cell_centered;
     amrex::ignore_unused(yp);
     i = static_cast<int>(amrex::Math::floor(x));
     W[0][1] = x - i;
 #   elif (defined WARPX_DIM_RZ)
-    const amrex::Real r = (std::sqrt(xp*xp+yp*yp) - plo[0]) * dxi[0] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
+    const amrex::Real r = (std::sqrt(xp*xp+yp*yp) - plo[0]) * dxi[0] - half_if_cell_centered;
     i = static_cast<int>(amrex::Math::floor(r));
     W[0][1] = r - i;
 #   endif
 
-    const amrex::Real z = (zp - plo[1]) * dxi[1] + static_cast<amrex::Real>(nodal-1)*0.5_rt;
+    const amrex::Real z = (zp - plo[1]) * dxi[1] - half_if_cell_centered;
     j = static_cast<int>(amrex::Math::floor(z));
     W[1][1] = z - j;
 
@@ -83,8 +92,15 @@ void compute_weights (const amrex::ParticleReal xp,
 
     k = 0;
 #else
-    amrex::ignore_unused(xp, yp, zp, plo, dxi, i, j, k, W, nodal);
-    ABLASTR_ABORT_WITH_MESSAGE("Error: compute_weights not yet implemented in 1D");
+    const amrex::Real z = (zp - plo[0]) * dxi[0] - half_if_cell_centered;
+    amrex::ignore_unused(xp, yp);
+    i = static_cast<int>(amrex::Math::floor(z));
+
+    W[0][1] = z - i;
+    W[0][0] = 1.0_rt - W[0][1];
+
+    j = 0;
+    k = 0;
 #endif
 }
 
@@ -100,8 +116,8 @@ amrex::Real interp_field_nodal (int i, int j, int k,
                                 const amrex::Real W[AMREX_SPACEDIM][2],
                                 amrex::Array4<const amrex::Real> const& scalar_field) noexcept
 {
-#if (defined WARPX_DIM_3D)
     amrex::Real value = 0;
+#if (defined WARPX_DIM_3D)
     value += scalar_field(i,   j  , k  ) * W[0][0] * W[1][0] * W[2][0];
     value += scalar_field(i+1, j  , k  ) * W[0][1] * W[1][0] * W[2][0];
     value += scalar_field(i,   j+1, k  ) * W[0][0] * W[1][1] * W[2][0];
@@ -111,15 +127,13 @@ amrex::Real interp_field_nodal (int i, int j, int k,
     value += scalar_field(i  , j+1, k+1) * W[0][0] * W[1][1] * W[2][1];
     value += scalar_field(i+1, j+1, k+1) * W[0][1] * W[1][1] * W[2][1];
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-    amrex::Real value = 0;
     value += scalar_field(i,   j ,  k) * W[0][0] * W[1][0];
     value += scalar_field(i+1, j ,  k) * W[0][1] * W[1][0];
     value += scalar_field(i,   j+1, k) * W[0][0] * W[1][1];
     value += scalar_field(i+1, j+1, k) * W[0][1] * W[1][1];
 #else
-    const amrex::Real value = 0;
-    amrex::ignore_unused(i, j, k, W, scalar_field);
-    ABLASTR_ABORT_WITH_MESSAGE("Error: interp_field not yet implemented in 1D");
+    value += scalar_field(i,   j ,  k) * W[0][0];
+    value += scalar_field(i+1, j ,  k) * W[0][1];
 #endif
     return value;
 }
@@ -144,7 +158,7 @@ amrex::Real doGatherScalarFieldNodal (const amrex::ParticleReal xp,
     // first find the weight of surrounding nodes to use during interpolation
     int ii, jj, kk;
     amrex::Real W[AMREX_SPACEDIM][2];
-    compute_weights(xp, yp, zp, lo, dxi, ii, jj, kk, W);
+    compute_weights<amrex::IndexType::NODE>(xp, yp, zp, lo, dxi, ii, jj, kk, W);
 
     return interp_field_nodal(ii, jj, kk, W, scalar_field);
 }
@@ -172,7 +186,7 @@ doGatherVectorFieldNodal (const amrex::ParticleReal xp,
     // first find the weight of surrounding nodes to use during interpolation
     int ii, jj, kk;
     amrex::Real W[AMREX_SPACEDIM][2];
-    compute_weights(xp, yp, zp, lo, dxi, ii, jj, kk, W);
+    compute_weights<amrex::IndexType::NODE>(xp, yp, zp, lo, dxi, ii, jj, kk, W);
 
     amrex::GpuArray<amrex::Real, 3> const field_interp = {
         interp_field_nodal(ii, jj, kk, W, vector_field_x),