Remove persistent E+B (#1050)

* add functor for doing the tmp particles copy for the back-transformed diagnosti * merge the particle push options into one kernel * EOL * fix assertion * add a FieldGatherandPushPX method to PhysicalParticleContainer * handle offset in copyAttribs * allow this functor to be constructed even it we aren't doing the back transformed diagnostics * EOL * update the overloads of PushPX for the Photon and RigidInjected ParticleContainers * function for dispatching the right field gather * init this val to 0.0 * fix some typos * handle scaling the fields for rigid injection * EOL * don't need to get pointers to E and B arrays in PushPX any more. * actually I can't remove these yet * EOL * variable order bug * move the QED stuff to the proper place * EOL * make sure we don't build these functors unless the runtime options are toggled * EOL * perform the field gather prior to the photon particle push * remove E and B components and FieldGather methods. Reimplement PushP for rigid injected and physical particles * update ionization to do field gather inline * remove E and B from the particle diagnostics * don't write E or B in these tests for particles * add missing files * remove EB from the Regtest ini file too * no need to do this twice * important typo * also do the gather inline for the QED processes that need to * move these sources inside ifdef for QED * fix bug in RZ * remove some fields from the Python tests. * remove all particle E and B comps from json benchmarks * don't assert that Ey is the langmuir output * remove uy from this output * update test * restore the mesh fields I turned off by mistake * turn off field IO for a few python tests I missed * fix typo * reset Langmuir_multi benchmark * update Langmuir_multi_nodal benchmark * update single precision langmuir bench * update psatd single precision languir one too * also do ionization_lab * finally, ionizaiton_boost * update benchmarks_json/Langmuir_multi_psatd.json * update benchmarks_json/Langmuir_multi_psatd_current_correction.json * update benchmarks_json/Langmuir_multi_psatd_momentum_conserving.json * update benchmarks_json/Langmuir_multi_psatd_nodal.json * remove the particle E and B from the choices in the docs * fix offset bug * also add the Gather subdirectory * Update Source/WarpX.H Co-authored-by: MaxThevenet <mthevenet@lbl.gov> * add docstring for LowerCornerWithGalilean * add new source files to CMakeLists.txt * also need to update the GPU regression tests * update the name of the output file for this python test * remove field gather call from FieldDiagnostics * fix typo in docstring * init fields to 0 * add docstring to the CopyParticleAttribs constructor * some explicit amrex::namepace Co-authored-by: MaxThevenet <mthevenet@lbl.gov>
author: Andrew Myers <atmyers@lbl.gov> 2020-06-22 09:28:16 -0700
committer: GitHub <noreply@github.com> 2020-06-22 09:28:16 -0700
commit: f4978e1001494e2b148128380fb37f3c2450209f (patch)
tree: 237a671ca194338411f16f52797e70114b220384
parent: 97755d1c2e04e5b8c3295182eaff472c73cf8a53 (diff)
download: WarpX-f4978e1001494e2b148128380fb37f3c2450209f.tar.gz
WarpX-f4978e1001494e2b148128380fb37f3c2450209f.tar.zst
WarpX-f4978e1001494e2b148128380fb37f3c2450209f.zip
101 files changed, 1517 insertions, 1510 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst
index 214ad8fe9..ab1634efe 100644
--- a/Docs/source/running_cpp/parameters.rst
+++ b/Docs/source/running_cpp/parameters.rst
@@ -1281,10 +1281,6 @@ In-situ capabilities can be used by turning on Sensei or Ascent (provided they a
 
     * ``ux`` ``uy`` ``uz`` for the particle momentum,
 
-    * ``Ex`` ``Ey`` ``Ez`` for the electric field on particles,
-
-    * ``Bx`` ``By`` ``Bz`` for the magnetic field on particles.
-
     The particle positions are always included.
     Use ``<species>.variables = none`` to plot no particle data, except particle position.
 
diff --git a/Examples/Modules/gaussian_beam/PICMI_inputs_gaussian_beam.py b/Examples/Modules/gaussian_beam/PICMI_inputs_gaussian_beam.py
index 737c5e94e..1e8851dd4 100644
--- a/Examples/Modules/gaussian_beam/PICMI_inputs_gaussian_beam.py
+++ b/Examples/Modules/gaussian_beam/PICMI_inputs_gaussian_beam.py
@@ -53,7 +53,7 @@ field_diag1 = picmi.FieldDiagnostic(name = 'diag1',
 part_diag1 = picmi.ParticleDiagnostic(name = 'diag1',
                                       period = 10,
                                       species = [electrons, protons],
-                                      data_list = ['weighting', 'momentum', 'fields'])
+                                      data_list = ['weighting', 'momentum'])
 
 sim = picmi.Simulation(solver = solver,
                        max_steps = 10,
diff --git a/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration.py b/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration.py
index 7fad60f19..07662ce3a 100644
--- a/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration.py
+++ b/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration.py
@@ -61,7 +61,7 @@ field_diag = picmi.FieldDiagnostic(name = 'diag1',
 part_diag = picmi.ParticleDiagnostic(name = 'diag1',
                                      period = max_steps,
                                      species = [beam, plasma],
-                                     data_list = ['ux', 'uy', 'uz', 'weighting', 'Ex', 'Ey', 'Ez'])
+                                     data_list = ['ux', 'uy', 'uz', 'weighting'])
 
 sim.add_diagnostic(field_diag)
 sim.add_diagnostic(part_diag)
diff --git a/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_mr.py b/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_mr.py
index bb0377e1a..665c78bc7 100644
--- a/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_mr.py
+++ b/Examples/Physics_applications/plasma_acceleration/PICMI_inputs_plasma_acceleration_mr.py
@@ -67,7 +67,7 @@ field_diag = picmi.FieldDiagnostic(name = 'diag1',
 part_diag = picmi.ParticleDiagnostic(name = 'diag1',
                                      period = 2,
                                      species = [beam, plasma],
-                                     data_list = ['ux', 'uy', 'uz', 'weighting', 'Ex', 'Ey', 'Ez'])
+                                     data_list = ['ux', 'uy', 'uz', 'weighting'])
 
 sim.add_diagnostic(field_diag)
 sim.add_diagnostic(part_diag)
diff --git a/Examples/Physics_applications/uniform_plasma/inputs_3d b/Examples/Physics_applications/uniform_plasma/inputs_3d
index 587e710f9..a8e3f14ec 100644
--- a/Examples/Physics_applications/uniform_plasma/inputs_3d
+++ b/Examples/Physics_applications/uniform_plasma/inputs_3d
@@ -38,7 +38,7 @@ electrons.uz_th  = 0.01 # uth the std of the (unitless) momentum
 diagnostics.diags_names = diag1 chk
 diag1.period = 4
 diag1.diag_type = Full
-diag1.electrons.variables = ux uy uz Ey Ez w
+diag1.electrons.variables = ux uy uz w
 
 chk.period = 6
 chk.diag_type = Full
diff --git a/Examples/Tests/Langmuir/PICMI_inputs_langmuir2d.py b/Examples/Tests/Langmuir/PICMI_inputs_langmuir2d.py
index 9d3dcc1db..3b67e990d 100644
--- a/Examples/Tests/Langmuir/PICMI_inputs_langmuir2d.py
+++ b/Examples/Tests/Langmuir/PICMI_inputs_langmuir2d.py
@@ -70,7 +70,7 @@ field_diag1 = picmi.FieldDiagnostic(name = 'diag1',
 part_diag1 = picmi.ParticleDiagnostic(name = 'diag1',
                                       period = diagnostic_interval,
                                       species = [electrons],
-                                      data_list = ['weighting', 'ux', 'Ex'])
+                                      data_list = ['weighting', 'ux'])
 
 ##########################
 # simulation setup
diff --git a/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rt.py b/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rt.py
index 8396f06a5..e6fd321cd 100644
--- a/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rt.py
+++ b/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rt.py
@@ -72,7 +72,7 @@ field_diag1 = picmi.FieldDiagnostic(name = 'diag1',
 part_diag1 = picmi.ParticleDiagnostic(name = 'diag1',
                                       period = diagnostic_interval,
                                       species = [electrons],
-                                      data_list = ['weighting', 'ux', 'Ex'])
+                                      data_list = ['weighting', 'ux'])
 
 ##########################
 # simulation setup
diff --git a/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rz_multimode_analyze.py b/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rz_multimode_analyze.py
index 37aad7308..e820a87db 100644
--- a/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rz_multimode_analyze.py
+++ b/Examples/Tests/Langmuir/PICMI_inputs_langmuir_rz_multimode_analyze.py
@@ -105,7 +105,7 @@ field_diag1 = picmi.FieldDiagnostic(name = 'diag1',
 part_diag1 = picmi.ParticleDiagnostic(name = 'diag1',
                                       period = diagnostic_interval,
                                       species = [electrons],
-                                      data_list = ['weighting', 'momentum', 'fields'])
+                                      data_list = ['weighting', 'momentum'])
 
 ##########################
 # simulation setup
diff --git a/Examples/Tests/Langmuir/analysis_langmuir_multi.py b/Examples/Tests/Langmuir/analysis_langmuir_multi.py
index f06aa4877..05fc99812 100755
--- a/Examples/Tests/Langmuir/analysis_langmuir_multi.py
+++ b/Examples/Tests/Langmuir/analysis_langmuir_multi.py
@@ -88,11 +88,8 @@ for field in ['particle_momentum_y',
     print('assert that this is NOT in ds.field_list', (species, field))
     assert (species, field) not in ds.field_list
 species = 'positrons'
-for field in ['particle_Ey']:
-    print('assert that this is in ds.field_list', (species, field))
-    assert (species, field) in ds.field_list
-for field in ['particle_momentum_y',
-              'particle_momentum_z']:
+for field in ['particle_momentum_x',
+              'particle_momentum_y']:
     print('assert that this is NOT in ds.field_list', (species, field))
     assert (species, field) not in ds.field_list
 
diff --git a/Examples/Tests/Langmuir/inputs_3d_multi_rt b/Examples/Tests/Langmuir/inputs_3d_multi_rt
index c9ddad495..708ee2a3e 100644
--- a/Examples/Tests/Langmuir/inputs_3d_multi_rt
+++ b/Examples/Tests/Langmuir/inputs_3d_multi_rt
@@ -88,4 +88,4 @@ diag1.period = 40
 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 = Ey
+diag1.positrons.variables = uz
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_2d.json b/Regression/Checksum/benchmarks_json/Langmuir_2d.json
index 85feaf801..ae115ffbe 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_2d.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_2d.json
@@ -1,6 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 22673325741654.633,
     "particle_cpu": 4096.0,
     "particle_id": 17567744.0,
     "particle_momentum_x": 9.644382989451675e-21,
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_2d_nompi.json b/Regression/Checksum/benchmarks_json/Langmuir_2d_nompi.json
index 921c28ccb..2d81edc17 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_2d_nompi.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_2d_nompi.json
@@ -1,6 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 22673325741654.633,
     "particle_cpu": 0.0,
     "particle_id": 52170752.0,
     "particle_momentum_x": 9.644382989451675e-21,
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_2d_single_precision.json b/Regression/Checksum/benchmarks_json/Langmuir_2d_single_precision.json
index 550106696..8e08ca0ed 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_2d_single_precision.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_2d_single_precision.json
@@ -1,6 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 22675062646016.0,
     "particle_cpu": 4096.0,
     "particle_id": 17305600.0,
     "particle_momentum_x": 9.646555269782835e-21,
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi.json b/Regression/Checksum/benchmarks_json/Langmuir_multi.json
index 3a0ef1456..bcb4b799b 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi.json
@@ -2,29 +2,29 @@
   "electrons": {
     "particle_cpu": 131072.0,
     "particle_id": 18862440448.0,
-    "particle_momentum_x": 9.638052142962567e-20,
+    "particle_momentum_x": 9.638052142962566e-20,
     "particle_position_x": 2.6214400000000015,
     "particle_position_y": 2.621440000000001,
     "particle_position_z": 2.6214400000000007,
     "particle_weight": 128000000000.00002
   },
   "lev=0": {
-    "Bx": 12.117994152440259,
-    "By": 12.117994153623144,
-    "Bz": 12.117994153640254,
-    "Ex": 84779179148604.14,
+    "Bx": 12.117994152442217,
+    "By": 12.117994153638133,
+    "Bz": 12.117994153639632,
+    "Ex": 84779179148604.16,
     "Ey": 84779179148604.05,
     "Ez": 84779179148604.05,
-    "jx": 6.087467475688615e+16,
-    "jy": 6.0874674756883144e+16,
-    "jz": 6.087467475688317e+16,
+    "jx": 6.087467475688619e+16,
+    "jy": 6.087467475688316e+16,
+    "jz": 6.087467475688315e+16,
     "part_per_cell": 524288.0,
-    "rho": 771603348.7307231
+    "rho": 771603348.7307227
   },
   "positrons": {
-    "particle_Ey": 84768126738155.34,
     "particle_cpu": 131072.0,
     "particle_id": 56518901760.0,
+    "particle_momentum_z": 9.638052142962866e-20,
     "particle_position_x": 2.6214400000000015,
     "particle_position_y": 2.621440000000001,
     "particle_position_z": 2.6214400000000007
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_nodal.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_nodal.json
index 346f5ed33..d3927b2ab 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_nodal.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_nodal.json
@@ -1,8 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 14986451875836.121,
-    "particle_Ey": 0.0,
-    "particle_Ez": 14986451875836.125,
     "particle_cpu": 32768.0,
     "particle_id": 1123057664.0,
     "particle_momentum_x": 5.668407425899705e-20,
@@ -24,9 +21,6 @@
     "jz": 1.0088447956385268e+16
   },
   "positrons": {
-    "particle_Ex": 14986451875836.117,
-    "particle_Ey": 0.0,
-    "particle_Ez": 14986451875836.125,
     "particle_cpu": 32768.0,
     "particle_id": 3371204608.0,
     "particle_momentum_x": 5.668407425899705e-20,
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd.json
index c3570d380..886342f01 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd.json
@@ -1,8 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 15085626458339.264,
-    "particle_Ey": 0.0,
-    "particle_Ez": 15085626458339.277,
     "particle_cpu": 32768.0,
     "particle_id": 1123057664.0,
     "particle_momentum_x": 5.665639562608442e-20,
@@ -22,9 +19,6 @@
     "part_per_cell": 131072.0
   },
   "positrons": {
-    "particle_Ex": 15085626458339.295,
-    "particle_Ey": 0.0,
-    "particle_Ez": 15085626458339.305,
     "particle_cpu": 32768.0,
     "particle_id": 3371204608.0,
     "particle_momentum_x": 5.665639562608429e-20,
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_current_correction.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_current_correction.json
index 52b54c8a3..34f48c585 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_current_correction.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_current_correction.json
@@ -1,8 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 15186050664939.146,
-    "particle_Ey": 0.0,
-    "particle_Ez": 15186050664939.232,
     "particle_cpu": 0.0,
     "particle_id": 2229305344.0,
     "particle_momentum_x": 5.658197625002215e-20,
@@ -24,9 +21,6 @@
     "rho": 21102031.294517133
   },
   "positrons": {
-    "particle_Ex": 15186050664939.146,
-    "particle_Ey": 0.0,
-    "particle_Ez": 15186050664939.232,
     "particle_cpu": 0.0,
     "particle_id": 6725599232.0,
     "particle_momentum_x": 5.658197625002216e-20,
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_momentum_conserving.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_momentum_conserving.json
index 662bb2366..063519515 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_momentum_conserving.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_momentum_conserving.json
@@ -1,8 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 14818473436129.969,
-    "particle_Ey": 0.0,
-    "particle_Ez": 14818473436129.988,
     "particle_cpu": 32768.0,
     "particle_id": 1123057664.0,
     "particle_momentum_x": 5.67711146356253e-20,
@@ -22,9 +19,6 @@
     "part_per_cell": 131072.0
   },
   "positrons": {
-    "particle_Ex": 14818473436129.98,
-    "particle_Ey": 0.0,
-    "particle_Ez": 14818473436130.002,
     "particle_cpu": 32768.0,
     "particle_id": 3371204608.0,
     "particle_momentum_x": 5.677111463562513e-20,
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_nodal.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_nodal.json
index 4c39d0e2c..ebef389de 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_nodal.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_nodal.json
@@ -1,8 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 14987240735672.16,
-    "particle_Ey": 0.0,
-    "particle_Ez": 14987238943607.992,
     "particle_cpu": 32768.0,
     "particle_id": 1123057664.0,
     "particle_momentum_x": 5.668629836343162e-20,
@@ -22,9 +19,6 @@
     "part_per_cell": 131072.0
   },
   "positrons": {
-    "particle_Ex": 14987240722150.396,
-    "particle_Ey": 0.0,
-    "particle_Ez": 14987238943288.441,
     "particle_cpu": 32768.0,
     "particle_id": 3371204608.0,
     "particle_momentum_x": 5.668629830474806e-20,
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_nodal.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_nodal.json
index 042d0c9fb..490566877 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_nodal.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_nodal.json
@@ -9,9 +9,9 @@
     "particle_weight": 128000000000.00002
   },
   "lev=0": {
-    "Bx": 17.67689485265776,
-    "By": 17.676894852668948,
-    "Bz": 17.676894852669562,
+    "Bx": 17.67689485265927,
+    "By": 17.676894852670383,
+    "Bz": 17.67689485267166,
     "Ex": 86079763548288.75,
     "Ey": 86079763548288.78,
     "Ez": 86079763548288.78,
@@ -19,12 +19,12 @@
     "jy": 5.803381905407015e+16,
     "jz": 5.803381905407016e+16,
     "part_per_cell": 524288.0,
-    "rho": 786759127.7574509
+    "rho": 786759127.7574512
   },
   "positrons": {
-    "particle_Ey": 86072707299548.95,
     "particle_cpu": 131072.0,
     "particle_id": 56518901760.0,
+    "particle_momentum_z": 9.320505028112306e-20,
     "particle_position_x": 2.6214400000000015,
     "particle_position_y": 2.621440000000001,
     "particle_position_z": 2.6214399999999998
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd.json
index e660cb6ee..7ef86212d 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd.json
@@ -2,29 +2,29 @@
   "electrons": {
     "particle_cpu": 131072.0,
     "particle_id": 18862440448.0,
-    "particle_momentum_x": 9.630000480732885e-20,
+    "particle_momentum_x": 9.630000480732875e-20,
     "particle_position_x": 2.621440000012802,
     "particle_position_y": 2.621440000012803,
     "particle_position_z": 2.6214399999999993,
     "particle_weight": 128000000000.00002
   },
   "lev=0": {
-    "Bx": 11.916179266713513,
-    "By": 11.91617926663706,
-    "Bz": 11.917607032204629,
-    "Ex": 84733730662264.23,
-    "Ey": 84733730662262.27,
-    "Ez": 84733724096318.2,
-    "jx": 6.0822131615524e+16,
-    "jy": 6.0822131615524344e+16,
-    "jz": 6.082209355386128e+16,
+    "Bx": 11.916179266720627,
+    "By": 11.916179266707902,
+    "Bz": 11.917607032225353,
+    "Ex": 84733730662264.19,
+    "Ey": 84733730662262.17,
+    "Ez": 84733724096318.14,
+    "jx": 6.082213161552392e+16,
+    "jy": 6.0822131615524264e+16,
+    "jz": 6.082209355386121e+16,
     "part_per_cell": 524288.0,
-    "rho": 702660848.6316689
+    "rho": 702660848.6316686
   },
   "positrons": {
-    "particle_Ey": 84722686943095.3,
     "particle_cpu": 131072.0,
     "particle_id": 56518901760.0,
+    "particle_momentum_z": 9.629994891851425e-20,
     "particle_position_x": 2.621440000012802,
     "particle_position_y": 2.6214400000128038,
     "particle_position_z": 2.6214399999999993
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_current_correction.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_current_correction.json
index cda05fce1..96141ff12 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_current_correction.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_current_correction.json
@@ -2,30 +2,30 @@
   "electrons": {
     "particle_cpu": 0.0,
     "particle_id": 37409652736.0,
-    "particle_momentum_x": 9.585443568545556e-20,
+    "particle_momentum_x": 9.585443568545559e-20,
     "particle_position_x": 2.6214400000000015,
     "particle_position_y": 2.621440000000001,
     "particle_position_z": 2.6214400000000007,
     "particle_weight": 128000000000.00002
   },
   "lev=0": {
-    "Bx": 11.867039048390716,
-    "By": 11.867039053074166,
-    "Bz": 11.867039053083186,
-    "Ex": 85001549508327.45,
-    "Ey": 85001549508324.98,
-    "Ez": 85001549508324.94,
+    "Bx": 11.867039048376148,
+    "By": 11.86703905304863,
+    "Bz": 11.867039053065115,
+    "Ex": 85001549508327.53,
+    "Ey": 85001549508324.97,
+    "Ez": 85001549508324.9,
     "divE": 7.97321190060971e+19,
-    "jx": 6.039975332058757e+16,
-    "jy": 6.0399753320588824e+16,
-    "jz": 6.039975332058885e+16,
+    "jx": 6.039975332058763e+16,
+    "jy": 6.039975332058887e+16,
+    "jz": 6.039975332058881e+16,
     "part_per_cell": 524288.0,
-    "rho": 705963156.392504
+    "rho": 705963156.3925041
   },
   "positrons": {
-    "particle_Ey": 84990471098283.88,
     "particle_cpu": 0.0,
     "particle_id": 112722575360.0,
+    "particle_momentum_z": 9.585443568545593e-20,
     "particle_position_x": 2.6214400000000015,
     "particle_position_y": 2.621440000000001,
     "particle_position_z": 2.6214400000000007
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_momentum_conserving.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_momentum_conserving.json
index dd5243ea2..4bd6b5e74 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_momentum_conserving.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_momentum_conserving.json
@@ -2,29 +2,29 @@
   "electrons": {
     "particle_cpu": 131072.0,
     "particle_id": 18862440448.0,
-    "particle_momentum_x": 9.313026326765913e-20,
+    "particle_momentum_x": 9.313026326765904e-20,
     "particle_position_x": 2.6214400000125964,
     "particle_position_y": 2.621440000012595,
     "particle_position_z": 2.62144,
     "particle_weight": 128000000000.00002
   },
   "lev=0": {
-    "Bx": 15.076849705026993,
-    "By": 15.076849706218196,
-    "Bz": 15.077115911503535,
-    "Ex": 86034772984177.27,
-    "Ey": 86034772984177.08,
-    "Ez": 86034760818989.12,
-    "jx": 5.854806321272413e+16,
-    "jy": 5.85480632127243e+16,
-    "jz": 5.854804280396043e+16,
+    "Bx": 15.076849705018862,
+    "By": 15.076849706227831,
+    "Bz": 15.07711591152293,
+    "Ex": 86034772984177.19,
+    "Ey": 86034772984177.0,
+    "Ez": 86034760818989.11,
+    "jx": 5.854806321272408e+16,
+    "jy": 5.854806321272421e+16,
+    "jz": 5.854804280396036e+16,
     "part_per_cell": 524288.0,
-    "rho": 720374346.6048291
+    "rho": 720374346.6048285
   },
   "positrons": {
-    "particle_Ey": 86027721802307.64,
     "particle_cpu": 131072.0,
     "particle_id": 56518901760.0,
+    "particle_momentum_z": 9.313023242196424e-20,
     "particle_position_x": 2.6214400000125964,
     "particle_position_y": 2.6214400000125946,
     "particle_position_z": 2.62144
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_nodal.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_nodal.json
index e394c6686..9230d0617 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_nodal.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_nodal.json
@@ -2,29 +2,29 @@
   "electrons": {
     "particle_cpu": 131072.0,
     "particle_id": 18862440448.0,
-    "particle_momentum_x": 9.336540089113062e-20,
+    "particle_momentum_x": 9.336540089113076e-20,
     "particle_position_x": 2.6214400000000015,
     "particle_position_y": 2.621440000000001,
     "particle_position_z": 2.621439999999999,
     "particle_weight": 128000000000.00002
   },
   "lev=0": {
-    "Bx": 17.35580092401342,
-    "By": 17.355800924024127,
-    "Bz": 17.35580092400305,
-    "Ex": 86173786000018.14,
-    "Ey": 86173786000018.19,
-    "Ez": 86173786000018.17,
-    "jx": 5.813715216795644e+16,
-    "jy": 5.813715216795643e+16,
-    "jz": 5.8137152167956424e+16,
+    "Bx": 17.355800924008413,
+    "By": 17.35580092399946,
+    "Bz": 17.355800923986315,
+    "Ex": 86173786000018.25,
+    "Ey": 86173786000018.25,
+    "Ez": 86173786000018.31,
+    "jx": 5.813715216795657e+16,
+    "jy": 5.813715216795654e+16,
+    "jz": 5.813715216795654e+16,
     "part_per_cell": 524288.0,
-    "rho": 721399745.4594865
+    "rho": 721399745.4594866
   },
   "positrons": {
-    "particle_Ey": 86166716262545.0,
     "particle_cpu": 131072.0,
     "particle_id": 56518901760.0,
+    "particle_momentum_z": 9.336540089113066e-20,
     "particle_position_x": 2.6214400000000015,
     "particle_position_y": 2.621440000000001,
     "particle_position_z": 2.621439999999999
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_single_precision.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_single_precision.json
index a0909be22..412d813f7 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_single_precision.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_psatd_single_precision.json
@@ -2,31 +2,31 @@
   "electrons": {
     "particle_cpu": 131072.0,
     "particle_id": 18270519296.0,
-    "particle_momentum_x": 9.630175835453285e-20,
-    "particle_position_x": 2.6214399486056834,
-    "particle_position_y": 2.6214399661659513,
-    "particle_position_z": 2.621439966153872,
+    "particle_momentum_x": 9.630167257760452e-20,
+    "particle_position_x": 2.6214399487543574,
+    "particle_position_y": 2.621439966081766,
+    "particle_position_z": 2.6214399662318897,
     "particle_weight": 127999991808.0
   },
   "lev=0": {
-    "Bx": 16.283409055462556,
-    "By": 27.847420256523446,
-    "Bz": 27.860402087754352,
-    "Ex": 84734483985576.88,
-    "Ey": 84733344104538.75,
-    "Ez": 84733290264251.75,
-    "jx": 6.08243835275767e+16,
-    "jy": 6.082200882039661e+16,
-    "jz": 6.082200213864928e+16,
+    "Bx": 16.292603092750593,
+    "By": 27.85642119481247,
+    "Bz": 27.88201210953058,
+    "Ex": 84734460310152.62,
+    "Ey": 84733293625145.25,
+    "Ez": 84733260535171.5,
+    "jx": 6.082428965504698e+16,
+    "jy": 6.082194972011667e+16,
+    "jz": 6.082192140821043e+16,
     "part_per_cell": 524288.0,
-    "rho": 702644169.0683594
+    "rho": 702644003.3183594
   },
   "positrons": {
-    "particle_Ey": 84722307724362.25,
     "particle_cpu": 131072.0,
     "particle_id": 54811295744.0,
-    "particle_position_x": 2.6214399497206955,
-    "particle_position_y": 2.6214399663009544,
-    "particle_position_z": 2.6214399657276033
+    "particle_momentum_z": 9.630067927781725e-20,
+    "particle_position_x": 2.6214399495141265,
+    "particle_position_y": 2.621439966217963,
+    "particle_position_z": 2.6214399657132788
   }
 }
 \ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_rz.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_rz.json
index e5e311f83..a2d7872d8 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_rz.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_rz.json
@@ -1,10 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 58.8416263852574,
-    "particle_By": 59.16597702799425,
-    "particle_Ex": 6600700367966.209,
-    "particle_Ey": 6543561338956.482,
-    "particle_Ez": 14386276215464.854,
     "particle_cpu": 29440.0,
     "particle_id": 1006524160.0,
     "particle_momentum_x": 3.3268808500279122e-21,
@@ -16,11 +11,6 @@
     "particle_weight": 81147583679.15044
   },
   "ions": {
-    "particle_Bx": 59.15077093467293,
-    "particle_By": 58.83567443208132,
-    "particle_Ex": 6556437015881.182,
-    "particle_Ey": 6595024813707.55,
-    "particle_Ez": 14386446938764.07,
     "particle_cpu": 29440.0,
     "particle_id": 3026343680.0,
     "particle_momentum_x": 1.6904985179025105e-20,
diff --git a/Regression/Checksum/benchmarks_json/Langmuir_multi_single_precision.json b/Regression/Checksum/benchmarks_json/Langmuir_multi_single_precision.json
index 1644dfadc..4b25682c9 100644
--- a/Regression/Checksum/benchmarks_json/Langmuir_multi_single_precision.json
+++ b/Regression/Checksum/benchmarks_json/Langmuir_multi_single_precision.json
@@ -2,31 +2,31 @@
   "electrons": {
     "particle_cpu": 131072.0,
     "particle_id": 18270519296.0,
-    "particle_momentum_x": 9.638160416482527e-20,
-    "particle_position_x": 2.6214399456764,
-    "particle_position_y": 2.6214399582924273,
-    "particle_position_z": 2.6214399582613623,
+    "particle_momentum_x": 9.638160826190717e-20,
+    "particle_position_x": 2.6214399456919466,
+    "particle_position_y": 2.6214399582432577,
+    "particle_position_z": 2.6214399583116688,
     "particle_weight": 127999991808.0
   },
   "lev=0": {
-    "Bx": 13.749213500551377,
-    "By": 23.70175167118765,
-    "Bz": 23.697557306031285,
-    "Ex": 84779697208316.75,
-    "Ey": 84779379998572.12,
-    "Ez": 84779381031286.25,
-    "jx": 6.087465586195686e+16,
-    "jy": 6.087430695259379e+16,
-    "jz": 6.087430637996634e+16,
+    "Bx": 13.746930533642427,
+    "By": 23.700993117762437,
+    "Bz": 23.69751037016036,
+    "Ex": 84779699984212.88,
+    "Ey": 84779380102544.12,
+    "Ez": 84779381453679.5,
+    "jx": 6.08746566247641e+16,
+    "jy": 6.087430753764646e+16,
+    "jz": 6.087430519585869e+16,
     "part_per_cell": 524288.0,
-    "rho": 771605195.3974609
+    "rho": 771605196.0458984
   },
   "positrons": {
-    "particle_Ey": 84768333819952.0,
     "particle_cpu": 131072.0,
     "particle_id": 54811295744.0,
-    "particle_position_x": 2.6214399459832123,
-    "particle_position_y": 2.621439958047347,
-    "particle_position_z": 2.6214399581733687
+    "particle_momentum_z": 9.638010595476651e-20,
+    "particle_position_x": 2.6214399459793185,
+    "particle_position_y": 2.621439958043453,
+    "particle_position_z": 2.6214399581214423
   }
 }
 \ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/Larmor.json b/Regression/Checksum/benchmarks_json/Larmor.json
index 0bed8ec6e..8c19e7b43 100644
--- a/Regression/Checksum/benchmarks_json/Larmor.json
+++ b/Regression/Checksum/benchmarks_json/Larmor.json
@@ -1,11 +1,5 @@
 {
   "electron": {
-    "particle_Bx": 0.0,
-    "particle_By": 0.0007789471980429488,
-    "particle_Bz": 0.0,
-    "particle_Ex": 19669.903641764264,
-    "particle_Ey": 0.0,
-    "particle_Ez": 537.4023205327452,
     "particle_cpu": 0.0,
     "particle_id": 1.0,
     "particle_momentum_x": 1.2084126370460285e-22,
@@ -44,12 +38,6 @@
     "part_per_cell": 0.0
   },
   "positron": {
-    "particle_Bx": 0.0,
-    "particle_By": 0.0007789471980429488,
-    "particle_Bz": 0.0,
-    "particle_Ex": 19669.90364176427,
-    "particle_Ey": 0.0,
-    "particle_Ez": 537.4023205327467,
     "particle_cpu": 0.0,
     "particle_id": 2.0,
     "particle_momentum_x": 1.2084126370460285e-22,
diff --git a/Regression/Checksum/benchmarks_json/LaserAcceleration.json b/Regression/Checksum/benchmarks_json/LaserAcceleration.json
index 27dac4b02..60f7e2a01 100644
--- a/Regression/Checksum/benchmarks_json/LaserAcceleration.json
+++ b/Regression/Checksum/benchmarks_json/LaserAcceleration.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 1823612.763420117,
-    "particle_By": 1594.036464819075,
-    "particle_Bz": 47860.82957589165,
-    "particle_Ex": 3400199636487.5806,
-    "particle_Ey": 726835961770741.9,
-    "particle_Ez": 34114413425319.66,
     "particle_cpu": 69212.0,
     "particle_id": 10198484142.0,
     "particle_momentum_x": 1.7927320949196398e-20,
diff --git a/Regression/Checksum/benchmarks_json/LaserAccelerationBoost.json b/Regression/Checksum/benchmarks_json/LaserAccelerationBoost.json
index c0adf0f38..d5e0cb243 100644
--- a/Regression/Checksum/benchmarks_json/LaserAccelerationBoost.json
+++ b/Regression/Checksum/benchmarks_json/LaserAccelerationBoost.json
@@ -1,11 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 2.0543216636907036,
-    "particle_By": 14.900417906604332,
-    "particle_Bz": 0.22233327031927602,
-    "particle_Ex": 4517976547.436662,
-    "particle_Ey": 633564460.6160239,
-    "particle_Ez": 15405468837.772211,
     "particle_cpu": 0.0,
     "particle_id": 35060500.0,
     "particle_momentum_x": 4.3244537833987973e-19,
@@ -16,12 +10,6 @@
     "particle_weight": 62415090744.60765
   },
   "electrons": {
-    "particle_Bx": 530.8105746418198,
-    "particle_By": 0.00010417687909189867,
-    "particle_Bz": 35.72584920547117,
-    "particle_Ex": 23373.585397386178,
-    "particle_Ey": 164837019826.8382,
-    "particle_Ez": 750684.6954455632,
     "particle_cpu": 1140.0,
     "particle_id": 43482810.0,
     "particle_momentum_x": 4.8513146378506466e-27,
@@ -32,12 +20,6 @@
     "particle_weight": 6.28658221073865e+16
   },
   "ions": {
-    "particle_Bx": 530.8105522074252,
-    "particle_By": 0.00010417687513907945,
-    "particle_Bz": 35.72584997453756,
-    "particle_Ex": 23373.584547084207,
-    "particle_Ey": 164837012558.85843,
-    "particle_Ez": 750684.6439845883,
     "particle_cpu": 1140.0,
     "particle_id": 43628730.0,
     "particle_momentum_x": 3.360051969962725e-30,
diff --git a/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json b/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json
index a19f55923..0e913059d 100644
--- a/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json
+++ b/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json
@@ -1,11 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 66.36298890249402,
-    "particle_By": 6824.872059894314,
-    "particle_Bz": 6.1501945107515,
-    "particle_Ex": 2030566419114.0784,
-    "particle_Ey": 21378773475.605545,
-    "particle_Ez": 1446901337250.283,
     "particle_cpu": 0.0,
     "particle_id": 214748364400.0,
     "particle_momentum_x": 4.3785337042075496e-20,
@@ -16,12 +10,6 @@
     "particle_weight": 12483018148921.525
   },
   "electrons": {
-    "particle_Bx": 8419020.374016589,
-    "particle_By": 1492.5583084263199,
-    "particle_Bz": 263998.6119992592,
-    "particle_Ex": 536859063153.5126,
-    "particle_Ey": 2485904866731919.0,
-    "particle_Ez": 2341935724590.825,
     "particle_cpu": 4788.0,
     "particle_id": 87752154.0,
     "particle_momentum_x": 4.275117840757358e-20,
diff --git a/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json b/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json
index 6ab7bd438..91a5cf2b8 100644
--- a/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json
+++ b/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json
@@ -1,10 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 1068.7120965886738,
-    "particle_By": 902.9886284146921,
-    "particle_Ex": 142251245678.64517,
-    "particle_Ey": 172725505795.35986,
-    "particle_Ez": 855052975051.5054,
     "particle_cpu": 0.0,
     "particle_id": 6865850.0,
     "particle_momentum_x": 4.3910959229136746e-20,
@@ -16,11 +11,6 @@
     "particle_weight": 6241509.074460764
   },
   "electrons": {
-    "particle_Bx": 33032.2917421148,
-    "particle_By": 179.35906104395352,
-    "particle_Ex": 60899865362.08995,
-    "particle_Ey": 9897446784083.63,
-    "particle_Ez": 205437656917.69284,
     "particle_cpu": 4128.0,
     "particle_id": 147732176.0,
     "particle_momentum_x": 3.2048286176059313e-24,
diff --git a/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost2d.json b/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost2d.json
index 8bbbf1631..0a4db7361 100644
--- a/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost2d.json
+++ b/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost2d.json
@@ -1,11 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 1.2421837369762971,
-    "particle_By": 112.61474112317613,
-    "particle_Bz": 0.15372744675077343,
-    "particle_Ex": 31266390602.63783,
-    "particle_Ey": 425776260.82294077,
-    "particle_Ez": 341412999657.737,
     "particle_cpu": 0.0,
     "particle_id": 18780500.0,
     "particle_momentum_x": 4.275414307401172e-19,
@@ -16,12 +10,6 @@
     "particle_weight": 3744905444676.4585
   },
   "driver": {
-    "particle_Bx": 0.17162902577074082,
-    "particle_By": 2599.284470679273,
-    "particle_Bz": 0.057532769068605646,
-    "particle_Ex": 383225305850.66095,
-    "particle_Ey": 56852018.90139305,
-    "particle_Ez": 3125011859579.5933,
     "particle_cpu": 0.0,
     "particle_id": 500500.0,
     "particle_momentum_x": 4.3244537833987973e-19,
@@ -43,12 +31,6 @@
     "jz": 102420439518871.45
   },
   "plasma_e": {
-    "particle_Bx": 0.0005857990454885307,
-    "particle_By": 1.404488613959929,
-    "particle_Bz": 0.00034742934066197536,
-    "particle_Ex": 406103675.3903916,
-    "particle_Ey": 201138.6662916025,
-    "particle_Ez": 204395172.69021654,
     "particle_cpu": 504.0,
     "particle_id": 8147196.0,
     "particle_momentum_x": 3.372310398663099e-24,
@@ -59,12 +41,6 @@
     "particle_weight": 677858073779801.2
   },
   "plasma_p": {
-    "particle_Bx": 0.0005858021375000984,
-    "particle_By": 1.4044879028870387,
-    "particle_Bz": 0.0003474290686070161,
-    "particle_Ex": 406103585.3352183,
-    "particle_Ey": 201139.68518534108,
-    "particle_Ez": 204396090.73847008,
     "particle_cpu": 504.0,
     "particle_id": 8806140.0,
     "particle_momentum_x": 3.3723091945226086e-24,
diff --git a/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost3d.json b/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost3d.json
index 00420b4f3..4144ea3df 100644
--- a/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost3d.json
+++ b/Regression/Checksum/benchmarks_json/PlasmaAccelerationBoost3d.json
@@ -1,11 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 3227.5589727359957,
-    "particle_By": 3088.1685105461247,
-    "particle_Bz": 4.57828126790789,
-    "particle_Ex": 309125717582.3938,
-    "particle_Ey": 328821757616.5381,
-    "particle_Ez": 6694077359865.185,
     "particle_cpu": 0.0,
     "particle_id": 550364500.0,
     "particle_momentum_x": 4.275414307401172e-19,
@@ -17,12 +11,6 @@
     "particle_weight": 3120754537.230381
   },
   "driver": {
-    "particle_Bx": 11659.346120355556,
-    "particle_By": 11215.848986117628,
-    "particle_Bz": 0.3070321484989495,
-    "particle_Ex": 530762394556.9341,
-    "particle_Ey": 552097091951.6976,
-    "particle_Ez": 5610847367213.652,
     "particle_cpu": 0.0,
     "particle_id": 500500.0,
     "particle_momentum_x": 4.747251762692475e+21,
@@ -34,12 +22,6 @@
     "particle_weight": 6241509074.460762
   },
   "driverback": {
-    "particle_Bx": 11688.980573229921,
-    "particle_By": 11159.188912235142,
-    "particle_Bz": 0.2544570753786498,
-    "particle_Ex": 523831138849.1566,
-    "particle_Ey": 551237711480.2166,
-    "particle_Ez": 5610277438559.053,
     "particle_cpu": 0.0,
     "particle_id": 551364500.0,
     "particle_momentum_x": 4.6831534563358655e+21,
@@ -62,12 +44,6 @@
     "jz": 1808527844770048.0
   },
   "plasma_e": {
-    "particle_Bx": 0.4955999326255563,
-    "particle_By": 0.5039645530718819,
-    "particle_Bz": 0.0001248192975851117,
-    "particle_Ex": 134686587.21869737,
-    "particle_Ey": 131645153.54607514,
-    "particle_Ez": 122964023.06369597,
     "particle_cpu": 3600.0,
     "particle_id": 1461156360.0,
     "particle_momentum_x": 2.0860691311063294e-24,
@@ -79,12 +55,6 @@
     "particle_weight": 3869103021.3870597
   },
   "plasma_p": {
-    "particle_Bx": 0.49560143627080033,
-    "particle_By": 0.5039657942943742,
-    "particle_Bz": 0.00012481965039337357,
-    "particle_Ex": 134687136.7309273,
-    "particle_Ey": 131645780.12750848,
-    "particle_Ez": 122964266.16215158,
     "particle_cpu": 3600.0,
     "particle_id": 1970570760.0,
     "particle_momentum_x": 2.0860707893867245e-24,
diff --git a/Regression/Checksum/benchmarks_json/PlasmaAccelerationMR.json b/Regression/Checksum/benchmarks_json/PlasmaAccelerationMR.json
index 145053393..955634bcb 100644
--- a/Regression/Checksum/benchmarks_json/PlasmaAccelerationMR.json
+++ b/Regression/Checksum/benchmarks_json/PlasmaAccelerationMR.json
@@ -1,11 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 32.31407180126963,
-    "particle_By": 6411.5062725212565,
-    "particle_Bz": 0.5571695868179536,
-    "particle_Ex": 1636006861658.8643,
-    "particle_Ey": 10062002545.305954,
-    "particle_Ez": 1935862949965.4546,
     "particle_cpu": 0.0,
     "particle_id": 2147483644000.0,
     "particle_momentum_x": 4.322181912633372e-19,
@@ -16,12 +10,6 @@
     "particle_weight": 12483018148921.53
   },
   "driver": {
-    "particle_Bx": 32.609326565006015,
-    "particle_By": 20031.368639484455,
-    "particle_Bz": 1.5411503609728023,
-    "particle_Ex": 5626282517515.15,
-    "particle_Ey": 9877108921.802748,
-    "particle_Ez": 18226789365193.43,
     "particle_cpu": 0.0,
     "particle_id": 10236500.0,
     "particle_momentum_x": 4.2756896805319457e-19,
@@ -54,12 +42,6 @@
     "jz": 4308416589615174.5
   },
   "plasma_e": {
-    "particle_Bx": 7.054183342141237,
-    "particle_By": 69097.79655524061,
-    "particle_Bz": 2.404677024126321,
-    "particle_Ex": 16213699604075.434,
-    "particle_Ey": 2276549918.021788,
-    "particle_Ez": 17965021901547.03,
     "particle_cpu": 3600.0,
     "particle_id": 41596488.0,
     "particle_momentum_x": 1.6052131138342657e-19,
diff --git a/Regression/Checksum/benchmarks_json/PlasmaMirror.json b/Regression/Checksum/benchmarks_json/PlasmaMirror.json
index bc59de5d0..620e5ad25 100644
--- a/Regression/Checksum/benchmarks_json/PlasmaMirror.json
+++ b/Regression/Checksum/benchmarks_json/PlasmaMirror.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 0.0,
-    "particle_By": 104961246.92730346,
-    "particle_Bz": 0.0,
-    "particle_Ex": 4.195946890380097e+16,
-    "particle_Ey": 0.0,
-    "particle_Ez": 5.398038858575869e+16,
     "particle_cpu": 3809.0,
     "particle_id": 115986415.0,
     "particle_momentum_x": 2.0003633196988347e-17,
@@ -16,12 +10,6 @@
     "particle_weight": 2.904173860599889e+18
   },
   "ions": {
-    "particle_Bx": 0.0,
-    "particle_By": 108988922.02755898,
-    "particle_Bz": 0.0,
-    "particle_Ex": 3.9621727637401224e+16,
-    "particle_Ey": 0.0,
-    "particle_Ez": 5.350781700255595e+16,
     "particle_cpu": 3834.0,
     "particle_id": 631174650.0,
     "particle_momentum_x": 2.805687685002513e-17,
diff --git a/Regression/Checksum/benchmarks_json/Python_Langmuir.json b/Regression/Checksum/benchmarks_json/Python_Langmuir.json
index 16fa4c9ab..16e593ccc 100644
--- a/Regression/Checksum/benchmarks_json/Python_Langmuir.json
+++ b/Regression/Checksum/benchmarks_json/Python_Langmuir.json
@@ -1,6 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 2.380866721344282e+16,
     "particle_cpu": 524288.0,
     "particle_id": 447927025664.0,
     "particle_momentum_x": 1.8915228266922305e-17,
diff --git a/Regression/Checksum/benchmarks_json/Python_Langmuir_2d.json b/Regression/Checksum/benchmarks_json/Python_Langmuir_2d.json
index bcb47cc19..0ebb255d1 100644
--- a/Regression/Checksum/benchmarks_json/Python_Langmuir_2d.json
+++ b/Regression/Checksum/benchmarks_json/Python_Langmuir_2d.json
@@ -1,6 +1,5 @@
 {
   "electrons": {
-    "particle_Ex": 219449752302051.72,
     "particle_cpu": 4096.0,
     "particle_id": 17567744.0,
     "particle_momentum_x": 1.0449418710231497e-19,
diff --git a/Regression/Checksum/benchmarks_json/Python_Langmuir_rz_multimode.json b/Regression/Checksum/benchmarks_json/Python_Langmuir_rz_multimode.json
index 4ce2699ec..8d55820e7 100644
--- a/Regression/Checksum/benchmarks_json/Python_Langmuir_rz_multimode.json
+++ b/Regression/Checksum/benchmarks_json/Python_Langmuir_rz_multimode.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 47.223888719248514,
-    "particle_By": 81.29252077224889,
-    "particle_Bz": 44.825833191000626,
-    "particle_Ex": 32416882559812.16,
-    "particle_Ey": 36617948944833.2,
-    "particle_Ez": 47879188522350.67,
     "particle_cpu": 353280.0,
     "particle_id": 158538771840.0,
     "particle_momentum_x": 2.1038307987216882e-20,
diff --git a/Regression/Checksum/benchmarks_json/Python_PlasmaAcceleration.json b/Regression/Checksum/benchmarks_json/Python_PlasmaAcceleration.json
index 1165cc1d3..fc079de81 100644
--- a/Regression/Checksum/benchmarks_json/Python_PlasmaAcceleration.json
+++ b/Regression/Checksum/benchmarks_json/Python_PlasmaAcceleration.json
@@ -1,8 +1,5 @@
 {
   "beam": {
-    "particle_Ex": 5756179858433.416,
-    "particle_Ey": 5756179858433.41,
-    "particle_Ez": 16964931026248.906,
     "particle_cpu": 0.0,
     "particle_id": 302285376.0,
     "particle_momentum_x": 5.708397145757236e-20,
@@ -23,9 +20,6 @@
     "part_per_cell": 590976.0
   },
   "plasma": {
-    "particle_Ex": 0.0,
-    "particle_Ey": 0.0,
-    "particle_Ez": 0.0,
     "particle_cpu": 589824.0,
     "particle_id": 514063630336.0,
     "particle_momentum_x": 0.0,
diff --git a/Regression/Checksum/benchmarks_json/Python_PlasmaAccelerationMR.json b/Regression/Checksum/benchmarks_json/Python_PlasmaAccelerationMR.json
index 89e3cbdb1..d3425a3ea 100644
--- a/Regression/Checksum/benchmarks_json/Python_PlasmaAccelerationMR.json
+++ b/Regression/Checksum/benchmarks_json/Python_PlasmaAccelerationMR.json
@@ -1,8 +1,5 @@
 {
   "beam": {
-    "particle_Ex": 107128787405.89337,
-    "particle_Ey": 107128787405.88715,
-    "particle_Ez": 118963821845577.23,
     "particle_cpu": 0.0,
     "particle_id": 46454366187008.0,
     "particle_momentum_x": 3.8119259477932654e-21,
@@ -32,9 +29,6 @@
     "part_per_cell": 21632.0
   },
   "plasma": {
-    "particle_Ex": 0.0,
-    "particle_Ey": 0.0,
-    "particle_Ez": 0.0,
     "particle_cpu": 8388608.0,
     "particle_id": 1231405243432960.0,
     "particle_momentum_x": 0.0,
diff --git a/Regression/Checksum/benchmarks_json/Python_gaussian_beam.json b/Regression/Checksum/benchmarks_json/Python_gaussian_beam.json
index fc66d5077..6b868e8ce 100644
--- a/Regression/Checksum/benchmarks_json/Python_gaussian_beam.json
+++ b/Regression/Checksum/benchmarks_json/Python_gaussian_beam.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 0.0002284627371099781,
-    "particle_By": 0.00021467707180826862,
-    "particle_Bz": 0.0002102644327743607,
-    "particle_Ex": 16981277.687556587,
-    "particle_Ey": 17009505.007470798,
-    "particle_Ez": 16939221.994517952,
     "particle_cpu": 0.0,
     "particle_id": 536887296.0,
     "particle_momentum_x": 6.797039213991982e-20,
@@ -29,12 +23,6 @@
     "part_per_cell": 65536.0
   },
   "protons": {
-    "particle_Bx": 0.00022692446773522987,
-    "particle_By": 0.00021360016225691481,
-    "particle_Bz": 0.00021031904641947448,
-    "particle_Ex": 16527170.45904002,
-    "particle_Ey": 16538759.196916314,
-    "particle_Ez": 16456286.211242234,
     "particle_cpu": 0.0,
     "particle_id": 1610629120.0,
     "particle_momentum_x": 3.199547530708718e-21,
diff --git a/Regression/Checksum/benchmarks_json/RigidInjection_lab.json b/Regression/Checksum/benchmarks_json/RigidInjection_lab.json
index 7b9fe6d7a..ee32c737e 100644
--- a/Regression/Checksum/benchmarks_json/RigidInjection_lab.json
+++ b/Regression/Checksum/benchmarks_json/RigidInjection_lab.json
@@ -1,11 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 2.6566894511996866e-05,
-    "particle_By": 0.0008575397894793725,
-    "particle_Bz": 3.651332298602814e-06,
-    "particle_Ex": 257689.1138980855,
-    "particle_Ey": 8160.730985703307,
-    "particle_Ez": 51350.93704570909,
     "particle_cpu": 0.0,
     "particle_id": 2001000.0,
     "particle_momentum_x": 4.2999340454016905e-17,
diff --git a/Regression/Checksum/benchmarks_json/TwoParticle_electrostatic.json b/Regression/Checksum/benchmarks_json/TwoParticle_electrostatic.json
index f1bcdcc2c..330cbc100 100644
--- a/Regression/Checksum/benchmarks_json/TwoParticle_electrostatic.json
+++ b/Regression/Checksum/benchmarks_json/TwoParticle_electrostatic.json
@@ -1,11 +1,5 @@
 {
   "electron1": {
-    "particle_Bx": 2.5036214403545607e-33,
-    "particle_By": 3.6007341350514865e-33,
-    "particle_Bz": 1.0868573384860435e-33,
-    "particle_Ex": 3.4123501299430056e-08,
-    "particle_Ey": 3.412350127746802e-08,
-    "particle_Ez": 3.4123501245407426e-08,
     "particle_cpu": 0.0,
     "particle_id": 1.0,
     "particle_momentum_x": 3.346088201352191e-29,
@@ -17,12 +11,6 @@
     "particle_weight": 1.0
   },
   "electron2": {
-    "particle_Bx": 1.96192602420799e-33,
-    "particle_By": 1.6876606788059644e-33,
-    "particle_Bz": 3.644858398013685e-33,
-    "particle_Ex": 3.4123501222575865e-08,
-    "particle_Ey": 3.412350127595063e-08,
-    "particle_Ez": 3.412350126765692e-08,
     "particle_cpu": 0.0,
     "particle_id": 2.0,
     "particle_momentum_x": 3.346088199424244e-29,
diff --git a/Regression/Checksum/benchmarks_json/Uniform_2d.json b/Regression/Checksum/benchmarks_json/Uniform_2d.json
index 0099c107f..caa7b033e 100644
--- a/Regression/Checksum/benchmarks_json/Uniform_2d.json
+++ b/Regression/Checksum/benchmarks_json/Uniform_2d.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 25700.400626747654,
-    "particle_By": 38491.27943906533,
-    "particle_Bz": 25938.261862056585,
-    "particle_Ex": 44859167061465.125,
-    "particle_Ey": 6794395639950.69,
-    "particle_Ez": 44201725613366.95,
     "particle_cpu": 32768.0,
     "particle_id": 1123057664.0,
     "particle_momentum_x": 1.351015959339252e-19,
diff --git a/Regression/Checksum/benchmarks_json/bilinear_filter.json b/Regression/Checksum/benchmarks_json/bilinear_filter.json
index 581fbd252..409e63bb2 100644
--- a/Regression/Checksum/benchmarks_json/bilinear_filter.json
+++ b/Regression/Checksum/benchmarks_json/bilinear_filter.json
@@ -1,11 +1,5 @@
 {
   "electron": {
-    "particle_Bx": 0.0,
-    "particle_By": 1.9043507547185993e-19,
-    "particle_Bz": 0.0,
-    "particle_Ex": 1.235067153470628e-09,
-    "particle_Ey": 0.0,
-    "particle_Ez": 0.0,
     "particle_cpu": 0.0,
     "particle_id": 1.0,
     "particle_momentum_x": 0.027309245307378237,
diff --git a/Regression/Checksum/benchmarks_json/galilean_2d_psatd.json b/Regression/Checksum/benchmarks_json/galilean_2d_psatd.json
index e00b5a517..444a33d87 100644
--- a/Regression/Checksum/benchmarks_json/galilean_2d_psatd.json
+++ b/Regression/Checksum/benchmarks_json/galilean_2d_psatd.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 0.0,
-    "particle_By": 0.009299012937832037,
-    "particle_Bz": 0.0,
-    "particle_Ex": 2816198.341358315,
-    "particle_Ey": 0.0,
-    "particle_Ez": 226824.31985528197,
     "particle_cpu": 0.0,
     "particle_id": 2160099328.0,
     "particle_momentum_x": 1.5318380864485464e-21,
@@ -16,12 +10,6 @@
     "particle_weight": 1.6888332018290936e+18
   },
   "ions": {
-    "particle_Bx": 0.0,
-    "particle_By": 0.00929911174906337,
-    "particle_Bz": 0.0,
-    "particle_Ex": 2816228.27643041,
-    "particle_Ey": 0.0,
-    "particle_Ez": 226827.1246175955,
     "particle_cpu": 0.0,
     "particle_id": 6522175488.0,
     "particle_momentum_x": 2.6093881778048826e-18,
diff --git a/Regression/Checksum/benchmarks_json/galilean_3d_psatd.json b/Regression/Checksum/benchmarks_json/galilean_3d_psatd.json
index d44296e24..d422efd30 100644
--- a/Regression/Checksum/benchmarks_json/galilean_3d_psatd.json
+++ b/Regression/Checksum/benchmarks_json/galilean_3d_psatd.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 8.686865792841713e-05,
-    "particle_By": 0.0002939590693725135,
-    "particle_Bz": 3.452350267327678e-05,
-    "particle_Ex": 88208.54084575009,
-    "particle_Ey": 28448.58165631292,
-    "particle_Ez": 7238.699716136691,
     "particle_cpu": 0.0,
     "particle_id": 524800.0,
     "particle_momentum_x": 2.7514767750837127e-23,
@@ -16,12 +10,6 @@
     "particle_weight": 1.0555207511431835e+17
   },
   "ions": {
-    "particle_Bx": 8.68715702754361e-05,
-    "particle_By": 0.0002939977457230219,
-    "particle_Bz": 3.452290180102611e-05,
-    "particle_Ex": 88220.42522764044,
-    "particle_Ey": 28448.864182536563,
-    "particle_Ez": 7238.973235144771,
     "particle_cpu": 0.0,
     "particle_id": 1573376.0,
     "particle_momentum_x": 4.033320884535886e-20,
diff --git a/Regression/Checksum/benchmarks_json/initial_distribution.json b/Regression/Checksum/benchmarks_json/initial_distribution.json
index 246d38041..967408203 100644
--- a/Regression/Checksum/benchmarks_json/initial_distribution.json
+++ b/Regression/Checksum/benchmarks_json/initial_distribution.json
@@ -1,11 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 1.7649330608090273e-13,
-    "particle_By": 1.6201399701460411e-13,
-    "particle_Bz": 1.66810333771523e-13,
-    "particle_Ex": 788306.6243597023,
-    "particle_Ey": 807373.0427832225,
-    "particle_Ez": 849117.4848964723,
     "particle_cpu": 0.0,
     "particle_id": 869690834940.0,
     "particle_momentum_x": 1.0643551806277579e-16,
@@ -17,12 +11,6 @@
     "particle_weight": 0.059552651788017566
   },
   "gaussian": {
-    "particle_Bx": 1.7020260238260468e-13,
-    "particle_By": 1.7933856122499433e-13,
-    "particle_Bz": 1.709281128949862e-13,
-    "particle_Ex": 849360.3198012591,
-    "particle_Ey": 875082.5083100101,
-    "particle_Ez": 878926.7763904883,
     "particle_cpu": 0.0,
     "particle_id": 131072256000.0,
     "particle_momentum_x": 1.1148694854418848e-18,
@@ -45,12 +33,6 @@
     "jz": 135925010961.82391
   },
   "maxwell_boltzmann": {
-    "particle_Bx": 1.7042273054862176e-13,
-    "particle_By": 1.7958882714803126e-13,
-    "particle_Bz": 1.7106339244943236e-13,
-    "particle_Ex": 849630.7492019709,
-    "particle_Ey": 874259.2333479257,
-    "particle_Ez": 879811.2999356263,
     "particle_cpu": 0.0,
     "particle_id": 393216256000.0,
     "particle_momentum_x": 1.1174920939750226e-18,
@@ -62,12 +44,6 @@
     "particle_weight": 8e+21
   },
   "maxwell_juttner": {
-    "particle_Bx": 1.7022872090850356e-13,
-    "particle_By": 1.7950700813504624e-13,
-    "particle_Bz": 1.709136282800873e-13,
-    "particle_Ex": 849523.6100404036,
-    "particle_Ey": 875528.2217687981,
-    "particle_Ez": 879351.2082461334,
     "particle_cpu": 0.0,
     "particle_id": 655360256000.0,
     "particle_momentum_x": 2.2132580272612802e-16,
diff --git a/Regression/Checksum/benchmarks_json/ionization_boost.json b/Regression/Checksum/benchmarks_json/ionization_boost.json
index f5f59ae46..8a1212864 100644
--- a/Regression/Checksum/benchmarks_json/ionization_boost.json
+++ b/Regression/Checksum/benchmarks_json/ionization_boost.json
@@ -1,35 +1,23 @@
 {
   "electrons": {
-    "particle_Bx": 0.0,
-    "particle_By": 22263824.862206914,
-    "particle_Bz": 0.0,
-    "particle_Ex": 9291877034473320.0,
-    "particle_Ey": 0.0,
-    "particle_Ez": 5319.164340861062,
-    "particle_cpu": 20.0,
-    "particle_id": 4332886865.0,
-    "particle_momentum_x": 2.0857482219894723e-17,
+    "particle_cpu": 27.0,
+    "particle_id": 4323919088.0,
+    "particle_momentum_x": 2.1096313075249046e-17,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 1.762015026993497e-17,
-    "particle_position_x": 0.11036837993083386,
-    "particle_position_y": 1.7283871495895655,
-    "particle_weight": 3.071068957031251e-09
+    "particle_momentum_z": 1.7264669463007245e-17,
+    "particle_position_x": 0.11048873689229222,
+    "particle_position_y": 1.7098255333343246,
+    "particle_weight": 3.0709523304687515e-09
   },
   "ions": {
-    "particle_Bx": 0.0,
-    "particle_By": 659186.9143798536,
-    "particle_Bz": 0.0,
-    "particle_Ex": 386286390605633.4,
-    "particle_Ey": 0.0,
-    "particle_Ez": 1130.542128053148,
     "particle_cpu": 0.0,
-    "particle_id": 587622528.0,
-    "particle_ionization_level": 52665.0,
-    "particle_momentum_x": 3.441120480256377e-18,
+    "particle_id": 588620608.0,
+    "particle_ionization_level": 52663.0,
+    "particle_momentum_x": 3.631023150973392e-18,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 1.0432995298680571e-13,
-    "particle_position_x": 0.02143988984627425,
-    "particle_position_y": 0.4726766035525105,
+    "particle_momentum_z": 1.043299530138721e-13,
+    "particle_position_x": 0.02143992896238411,
+    "particle_position_y": 0.47267660340847956,
     "particle_weight": 5.000947000000002e-10
   },
   "lev=0": {
@@ -38,9 +26,9 @@
     "Bz": 0.0,
     "Ex": 5472946208912828.0,
     "Ey": 0.0,
-    "Ez": 829.4153550071132,
+    "Ez": 829.4153550071181,
     "jx": 12255073032842.39,
     "jy": 0.0,
-    "jz": 83416.00000110811
+    "jz": 83416.00000112198
   }
 }
 \ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/ionization_lab.json b/Regression/Checksum/benchmarks_json/ionization_lab.json
index 5f0a437fd..0113e23df 100644
--- a/Regression/Checksum/benchmarks_json/ionization_lab.json
+++ b/Regression/Checksum/benchmarks_json/ionization_lab.json
@@ -1,35 +1,23 @@
 {
   "electrons": {
-    "particle_Bx": 0.0,
-    "particle_By": 160958142.7106404,
-    "particle_Bz": 0.0,
-    "particle_Ex": 4.836891522729804e+16,
-    "particle_Ey": 0.0,
-    "particle_Ez": 6681.416407321183,
-    "particle_cpu": 17886.0,
-    "particle_id": 1440973250.0,
-    "particle_momentum_x": 4.345383730592745e-18,
+    "particle_cpu": 17826.0,
+    "particle_id": 1434117982.0,
+    "particle_momentum_x": 4.423317346074595e-18,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 2.6255247164666637e-18,
-    "particle_position_x": 0.11067439437835992,
-    "particle_position_y": 0.6395301140823729,
-    "particle_weight": 3.45546875e-10
+    "particle_momentum_z": 2.654287388285012e-18,
+    "particle_position_x": 0.11034345761508463,
+    "particle_position_y": 0.6415606590060071,
+    "particle_weight": 3.440703125e-10
   },
   "ions": {
-    "particle_Bx": 0.0,
-    "particle_By": 10930527.873315271,
-    "particle_Bz": 0.0,
-    "particle_Ex": 3270566810756972.0,
-    "particle_Ey": 0.0,
-    "particle_Ez": 1887.110708409788,
     "particle_cpu": 4864.0,
     "particle_id": 193800448.0,
-    "particle_ionization_level": 72810.0,
-    "particle_momentum_x": 1.73229610377399e-18,
+    "particle_ionization_level": 72773.0,
+    "particle_momentum_x": 1.7335485753555193e-18,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 3.509872456624198e-23,
-    "particle_position_x": 0.032000004532743265,
-    "particle_position_y": 0.12800000460884142,
+    "particle_momentum_z": 3.5110738038450676e-23,
+    "particle_position_x": 0.03199996960212278,
+    "particle_position_y": 0.12800000461301925,
     "particle_weight": 9.999999999999999e-11
   },
   "lev=0": {
@@ -38,9 +26,9 @@
     "Bz": 0.0,
     "Ex": 7878018647090106.0,
     "Ey": 0.0,
-    "Ez": 2740.9364504831315,
+    "Ez": 2740.936450483133,
     "jx": 1.2111228466502686e+16,
     "jy": 0.0,
-    "jz": 1.360852807841737e-07
+    "jz": 1.362596673906476e-07
   }
 }
 \ No newline at end of file
diff --git a/Regression/Checksum/benchmarks_json/momentum-conserving-gather.json b/Regression/Checksum/benchmarks_json/momentum-conserving-gather.json
index 3d8525b25..2e6967051 100644
--- a/Regression/Checksum/benchmarks_json/momentum-conserving-gather.json
+++ b/Regression/Checksum/benchmarks_json/momentum-conserving-gather.json
@@ -1,11 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 30.446171922907517,
-    "particle_By": 3870.917493104042,
-    "particle_Bz": 0.34755106738400265,
-    "particle_Ex": 1014015529098.6418,
-    "particle_Ey": 19826252289.36772,
-    "particle_Ez": 1908690024876.2688,
     "particle_cpu": 0.0,
     "particle_id": 2147483644000.0,
     "particle_momentum_x": 4.323192428952614e-19,
@@ -16,12 +10,6 @@
     "particle_weight": 12483018148921.53
   },
   "driver": {
-    "particle_Bx": 31.940096497888742,
-    "particle_By": 15604.732356054985,
-    "particle_Bz": 1.2252572333665939,
-    "particle_Ex": 4364681742392.5117,
-    "particle_Ey": 9873744504.0512,
-    "particle_Ez": 18055631247141.04,
     "particle_cpu": 0.0,
     "particle_id": 10236500.0,
     "particle_momentum_x": 4.273571972654775e-19,
@@ -54,12 +42,6 @@
     "jz": 4439666503460565.0
   },
   "plasma_e": {
-    "particle_Bx": 6.9776518155443075,
-    "particle_By": 65971.4138126626,
-    "particle_Bz": 2.2431702261360305,
-    "particle_Ex": 14972474416018.527,
-    "particle_Ey": 2292688361.624074,
-    "particle_Ez": 17838346022653.637,
     "particle_cpu": 3600.0,
     "particle_id": 41596488.0,
     "particle_momentum_x": 1.5130313200302403e-19,
diff --git a/Regression/Checksum/benchmarks_json/particle_pusher.json b/Regression/Checksum/benchmarks_json/particle_pusher.json
index dd2548827..7c6c52032 100644
--- a/Regression/Checksum/benchmarks_json/particle_pusher.json
+++ b/Regression/Checksum/benchmarks_json/particle_pusher.json
@@ -11,12 +11,6 @@
     "jz": 0.0
   },
   "positron": {
-    "particle_Bx": 0.0,
-    "particle_By": 0.0,
-    "particle_Bz": 1.0,
-    "particle_Ex": 299417482.9214179,
-    "particle_Ey": 0.0,
-    "particle_Ez": 0.0,
     "particle_cpu": 0.0,
     "particle_id": 1.0,
     "particle_momentum_x": 6.55085314065218e-05,
diff --git a/Regression/Checksum/benchmarks_json/qed_breit_wheeler_opt_depth_evolution.json b/Regression/Checksum/benchmarks_json/qed_breit_wheeler_opt_depth_evolution.json
index 0e2c5e60f..5553f1917 100644
--- a/Regression/Checksum/benchmarks_json/qed_breit_wheeler_opt_depth_evolution.json
+++ b/Regression/Checksum/benchmarks_json/qed_breit_wheeler_opt_depth_evolution.json
@@ -1,11 +1,5 @@
 {
   "ele_bw": {
-    "particle_Bx": 44488571428.57164,
-    "particle_By": 66732857142.85782,
-    "particle_Bz": 133465714285.71455,
-    "particle_Ex": 3.7889246226388705e+19,
-    "particle_Ey": 1.5155698490565198e+19,
-    "particle_Ez": 2.273354773582188e+19,
     "particle_cpu": 7868.0,
     "particle_id": 909726920.0,
     "particle_momentum_x": 2.5371833905509472e-14,
@@ -64,12 +58,6 @@
     "particle_position_z": 0.0016801197916666664
   },
   "pos_bw": {
-    "particle_Bx": 44488571428.57158,
-    "particle_By": 66732857142.85768,
-    "particle_Bz": 133465714285.71454,
-    "particle_Ex": 3.788924622638887e+19,
-    "particle_Ey": 1.5155698490565136e+19,
-    "particle_Ez": 2.2733547735821787e+19,
     "particle_cpu": 7868.0,
     "particle_id": 911804493.0,
     "particle_momentum_x": 2.5161516094347012e-14,
diff --git a/Regression/Checksum/benchmarks_json/qed_breit_wheeler_tau_init.json b/Regression/Checksum/benchmarks_json/qed_breit_wheeler_tau_init.json
index 553f043fb..2241e0fb5 100644
--- a/Regression/Checksum/benchmarks_json/qed_breit_wheeler_tau_init.json
+++ b/Regression/Checksum/benchmarks_json/qed_breit_wheeler_tau_init.json
@@ -11,12 +11,6 @@
     "jz": 0.0
   },
   "photons": {
-    "particle_Bx": 0.0,
-    "particle_By": 0.0,
-    "particle_Bz": 0.0,
-    "particle_Ex": 0.0,
-    "particle_Ey": 0.0,
-    "particle_Ez": 0.0,
     "particle_cpu": 28800.0,
     "particle_id": 993882240.0,
     "particle_momentum_x": 1.257725971841766e-15,
diff --git a/Regression/Checksum/benchmarks_json/qed_quantum_sync_tau_init.json b/Regression/Checksum/benchmarks_json/qed_quantum_sync_tau_init.json
index 26133ef69..6196b8aa6 100644
--- a/Regression/Checksum/benchmarks_json/qed_quantum_sync_tau_init.json
+++ b/Regression/Checksum/benchmarks_json/qed_quantum_sync_tau_init.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 0.23333485603512907,
-    "particle_By": 0.3306506464928457,
-    "particle_Bz": 0.23170804111682167,
-    "particle_Ex": 306898947.1037948,
-    "particle_Ey": 314678746.66859293,
-    "particle_Ez": 322524922.65796685,
     "particle_cpu": 28800.0,
     "particle_id": 993882240.0,
     "particle_momentum_x": 1.2577259718364536e-15,
@@ -28,12 +22,6 @@
     "jz": 665944376473.1846
   },
   "positrons": {
-    "particle_Bx": 0.23364278767344948,
-    "particle_By": 0.3313944153647703,
-    "particle_Bz": 0.23186653729970977,
-    "particle_Ex": 307442525.02915114,
-    "particle_Ey": 313908687.7726069,
-    "particle_Ez": 322478933.98023397,
     "particle_cpu": 28800.0,
     "particle_id": 2969792640.0,
     "particle_momentum_x": 1.2564867705042074e-15,
diff --git a/Regression/Checksum/benchmarks_json/qed_schwinger2.json b/Regression/Checksum/benchmarks_json/qed_schwinger2.json
index cd3ecf6c9..e2a524929 100644
--- a/Regression/Checksum/benchmarks_json/qed_schwinger2.json
+++ b/Regression/Checksum/benchmarks_json/qed_schwinger2.json
@@ -1,11 +1,5 @@
 {
   "ele_schwinger": {
-    "particle_Bx": 3439183579241.8213,
-    "particle_By": 1076561948990.9731,
-    "particle_Bz": 3760851107750.2676,
-    "particle_Ex": 1.6138262991542085e+26,
-    "particle_Ey": 1.0214830999058134e+25,
-    "particle_Ez": 3.5684391886872975e+25,
     "particle_cpu": 1024.0,
     "particle_id": 1573888.0,
     "particle_momentum_x": 1.253856079274937e-09,
@@ -28,12 +22,6 @@
     "jz": 3.214948145600786e+30
   },
   "pos_schwinger": {
-    "particle_Bx": 3439183579241.8213,
-    "particle_By": 1076561948990.9731,
-    "particle_Bz": 3760851107750.2676,
-    "particle_Ex": 1.6138262991542085e+26,
-    "particle_Ey": 1.0214830999058134e+25,
-    "particle_Ez": 3.5684391886872975e+25,
     "particle_cpu": 1024.0,
     "particle_id": 2622464.0,
     "particle_momentum_x": 1.2705200005945566e-09,
diff --git a/Regression/Checksum/benchmarks_json/qed_schwinger3.json b/Regression/Checksum/benchmarks_json/qed_schwinger3.json
index ea687e711..d2de759a4 100644
--- a/Regression/Checksum/benchmarks_json/qed_schwinger3.json
+++ b/Regression/Checksum/benchmarks_json/qed_schwinger3.json
@@ -1,11 +1,5 @@
 {
   "ele_schwinger": {
-    "particle_Bx": 1.348034558144351,
-    "particle_By": 0.0,
-    "particle_Bz": 0.33927654383456507,
-    "particle_Ex": 0.0,
-    "particle_Ey": 1.3909415199145668e+20,
-    "particle_Ez": 0.0,
     "particle_cpu": 631.0,
     "particle_id": 610043.0,
     "particle_momentum_x": 7.391180107793151e-28,
@@ -28,12 +22,6 @@
     "jz": 0.0
   },
   "pos_schwinger": {
-    "particle_Bx": 1.3398489034143053,
-    "particle_By": 0.0,
-    "particle_Bz": 0.3427233099177752,
-    "particle_Ex": 0.0,
-    "particle_Ey": 1.3909415199145551e+20,
-    "particle_Ez": 0.0,
     "particle_cpu": 631.0,
     "particle_id": 1017026.0,
     "particle_momentum_x": 7.723864257907636e-28,
diff --git a/Regression/Checksum/benchmarks_json/qed_schwinger4.json b/Regression/Checksum/benchmarks_json/qed_schwinger4.json
index db5ad389e..b2276d547 100644
--- a/Regression/Checksum/benchmarks_json/qed_schwinger4.json
+++ b/Regression/Checksum/benchmarks_json/qed_schwinger4.json
@@ -1,11 +1,5 @@
 {
   "ele_schwinger": {
-    "particle_Bx": 2.319332155934254,
-    "particle_By": 1707847966720000.0,
-    "particle_Bz": 0.6015864770059552,
-    "particle_Ex": 2139883430.888652,
-    "particle_Ey": 0.0,
-    "particle_Ez": 5.119999999998239e+23,
     "particle_cpu": 1024.0,
     "particle_id": 1573888.0,
     "particle_momentum_x": 6.569424596732217e-12,
@@ -28,12 +22,6 @@
     "jz": 1.5866782361655698e+16
   },
   "pos_schwinger": {
-    "particle_Bx": 2.3212042547187046,
-    "particle_By": 1707847966720000.0,
-    "particle_Bz": 0.6075072929142599,
-    "particle_Ex": 2146992296.4054935,
-    "particle_Ey": 0.0,
-    "particle_Ez": 5.119999999998239e+23,
     "particle_cpu": 1024.0,
     "particle_id": 2622464.0,
     "particle_momentum_x": 6.569424596732212e-12,
diff --git a/Regression/Checksum/benchmarks_json/radiation_reaction.json b/Regression/Checksum/benchmarks_json/radiation_reaction.json
index 9466924a2..4bf975174 100644
--- a/Regression/Checksum/benchmarks_json/radiation_reaction.json
+++ b/Regression/Checksum/benchmarks_json/radiation_reaction.json
@@ -1,11 +1,5 @@
 {
   "ele_para0": {
-    "particle_Bx": 917978.2333474257,
-    "particle_By": 1376967.350021139,
-    "particle_Bz": 2753934.700042278,
-    "particle_Ex": 0.003842707741931274,
-    "particle_Ey": 0.015625799690926225,
-    "particle_Ez": 0.023446517674477488,
     "particle_cpu": 0.0,
     "particle_id": 1.0,
     "particle_momentum_x": 7.802641516393779e-20,
@@ -17,12 +11,6 @@
     "particle_weight": 1e-08
   },
   "ele_perp0": {
-    "particle_Bx": 917978.2333474257,
-    "particle_By": 1376967.350021139,
-    "particle_Bz": 2753934.700042278,
-    "particle_Ex": 0.1250313976349998,
-    "particle_Ey": 0.0776652150231756,
-    "particle_Ez": 0.06187021356968523,
     "particle_cpu": 0.0,
     "particle_id": 2.0,
     "particle_momentum_x": 8.11389431708072e-21,
@@ -34,12 +22,6 @@
     "particle_weight": 1e-08
   },
   "ele_perp1": {
-    "particle_Bx": 917978.2333474257,
-    "particle_By": 1376967.350021139,
-    "particle_Bz": 2753934.700042278,
-    "particle_Ex": 0.03710062023180817,
-    "particle_Ey": 0.002435032603780963,
-    "particle_Ez": 0.023704514044336038,
     "particle_cpu": 0.0,
     "particle_id": 3.0,
     "particle_momentum_x": 2.3334382603801444e-20,
@@ -51,12 +33,6 @@
     "particle_weight": 1e-08
   },
   "ele_perp2": {
-    "particle_Bx": 917978.2333474257,
-    "particle_By": 1376967.350021139,
-    "particle_Bz": 2753934.700042278,
-    "particle_Ex": 0.0002716012778477919,
-    "particle_Ey": 0.005385504253127432,
-    "particle_Ez": 0.016067456901663155,
     "particle_cpu": 0.0,
     "particle_id": 4.0,
     "particle_momentum_x": 1.81896887291916e-20,
@@ -71,12 +47,6 @@
     "rho": 0.00051269652288
   },
   "pos_perp2": {
-    "particle_Bx": 917978.2333474257,
-    "particle_By": 1376967.350021139,
-    "particle_Bz": 2753934.700042278,
-    "particle_Ex": 0.004491380475289906,
-    "particle_Ey": 0.012777393943318885,
-    "particle_Ez": 0.02007316097507273,
     "particle_cpu": 0.0,
     "particle_id": 5.0,
     "particle_momentum_x": 1.818968872919164e-20,
diff --git a/Regression/Checksum/benchmarks_json/reduced_diags.json b/Regression/Checksum/benchmarks_json/reduced_diags.json
index 023995db0..8a385327d 100644
--- a/Regression/Checksum/benchmarks_json/reduced_diags.json
+++ b/Regression/Checksum/benchmarks_json/reduced_diags.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 0.08468252959119768,
-    "particle_By": 0.08433759166060013,
-    "particle_Bz": 0.0834983764785095,
-    "particle_Ex": 106389075.07051627,
-    "particle_Ey": 110696689.8493621,
-    "particle_Ez": 107474460.40560272,
     "particle_cpu": 16384.0,
     "particle_id": 268451840.0,
     "particle_momentum_x": 2.433513100910759e-19,
@@ -28,12 +22,6 @@
     "jz": 719566.2651192012
   },
   "photons": {
-    "particle_Bx": 0.08443156833972876,
-    "particle_By": 0.08498460548568723,
-    "particle_Bz": 0.08428210467483634,
-    "particle_Ex": 106054332.53718573,
-    "particle_Ey": 111005174.77433285,
-    "particle_Ez": 105476983.07484898,
     "particle_cpu": 16384.0,
     "particle_id": 1342193664.0,
     "particle_momentum_x": 1.4328434977988592e-18,
@@ -45,12 +33,6 @@
     "particle_weight": 800000000000000.0
   },
   "protons": {
-    "particle_Bx": 0.08405060899990621,
-    "particle_By": 0.08395398257418159,
-    "particle_Bz": 0.08318247160660851,
-    "particle_Ex": 102197583.14430666,
-    "particle_Ey": 106376868.79571314,
-    "particle_Ez": 102627734.35036321,
     "particle_cpu": 16384.0,
     "particle_id": 805322752.0,
     "particle_momentum_x": 1.4305311394194743e-19,
diff --git a/Regression/Checksum/benchmarks_json/reduced_diags_loadbalancecosts_heuristic.json b/Regression/Checksum/benchmarks_json/reduced_diags_loadbalancecosts_heuristic.json
index 36eb4d421..650543410 100644
--- a/Regression/Checksum/benchmarks_json/reduced_diags_loadbalancecosts_heuristic.json
+++ b/Regression/Checksum/benchmarks_json/reduced_diags_loadbalancecosts_heuristic.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 0.0,
-    "particle_By": 0.0,
-    "particle_Bz": 0.0,
-    "particle_Ex": 0.0,
-    "particle_Ey": 0.0,
-    "particle_Ez": 0.0,
     "particle_cpu": 0.0,
     "particle_id": 8590000128.0,
     "particle_momentum_x": 0.0,
diff --git a/Regression/Checksum/benchmarks_json/reduced_diags_loadbalancecosts_timers.json b/Regression/Checksum/benchmarks_json/reduced_diags_loadbalancecosts_timers.json
index 36eb4d421..650543410 100644
--- a/Regression/Checksum/benchmarks_json/reduced_diags_loadbalancecosts_timers.json
+++ b/Regression/Checksum/benchmarks_json/reduced_diags_loadbalancecosts_timers.json
@@ -1,11 +1,5 @@
 {
   "electrons": {
-    "particle_Bx": 0.0,
-    "particle_By": 0.0,
-    "particle_Bz": 0.0,
-    "particle_Ex": 0.0,
-    "particle_Ey": 0.0,
-    "particle_Ez": 0.0,
     "particle_cpu": 0.0,
     "particle_id": 8590000128.0,
     "particle_momentum_x": 0.0,
diff --git a/Regression/Checksum/benchmarks_json/subcyclingMR.json b/Regression/Checksum/benchmarks_json/subcyclingMR.json
index ac83e2c46..c5564282d 100644
--- a/Regression/Checksum/benchmarks_json/subcyclingMR.json
+++ b/Regression/Checksum/benchmarks_json/subcyclingMR.json
@@ -1,11 +1,5 @@
 {
   "beam": {
-    "particle_Bx": 0.0,
-    "particle_By": 1295546.8975611941,
-    "particle_Bz": 0.0,
-    "particle_Ex": 344765556152058.25,
-    "particle_Ey": 0.0,
-    "particle_Ez": 1193101601971968.0,
     "particle_cpu": 0.0,
     "particle_id": 21474836440000.0,
     "particle_momentum_x": 3.0267885065250823e-19,
@@ -16,12 +10,6 @@
     "particle_weight": 62415090744607.65
   },
   "driver": {
-    "particle_Bx": 0.0,
-    "particle_By": 7036528.904102431,
-    "particle_Bz": 0.0,
-    "particle_Ex": 2076986724754339.0,
-    "particle_Ey": 0.0,
-    "particle_Ez": 1520191782419694.0,
     "particle_cpu": 0.0,
     "particle_id": 50005000.0,
     "particle_momentum_x": 7.231072222092701e-19,
@@ -54,12 +42,6 @@
     "jz": 3.914141701645454e+17
   },
   "plasma_e": {
-    "particle_Bx": 0.0,
-    "particle_By": 3453725.2247103425,
-    "particle_Bz": 0.0,
-    "particle_Ex": 1427604131682412.5,
-    "particle_Ey": 0.0,
-    "particle_Ez": 1484821487586464.8,
     "particle_cpu": 0.0,
     "particle_id": 22976894251680.0,
     "particle_momentum_x": 1.1529420490113797e-18,
@@ -70,12 +52,6 @@
     "particle_weight": 5723510742187501.0
   },
   "plasma_p": {
-    "particle_Bx": 0.0,
-    "particle_By": 3766950.28595219,
-    "particle_Bz": 0.0,
-    "particle_Ex": 1592516445686014.5,
-    "particle_Ey": 0.0,
-    "particle_Ez": 1810017527838155.5,
     "particle_cpu": 0.0,
     "particle_id": 23097613251264.0,
     "particle_momentum_x": 1.5163562656560446e-18,
diff --git a/Regression/WarpX-GPU-tests.ini b/Regression/WarpX-GPU-tests.ini
index b2fc21831..97a955753 100644
--- a/Regression/WarpX-GPU-tests.ini
+++ b/Regression/WarpX-GPU-tests.ini
@@ -236,7 +236,7 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 particleTypes = electrons
-runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux Ex
+runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir2d.py
 analysisOutputImage = langmuir2d_analysis.png
 tolerance = 1e-12
@@ -244,7 +244,7 @@ tolerance = 1e-12
 [Langmuir_2d_single_precision]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_rt
-runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux Ex
+runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux
 dim = 2
 addToCompileString = USE_GPU=TRUE PRECISION=FLOAT USE_SINGLE_PRECISION_PARTICLES=TRUE
 restartTest = 0
@@ -274,7 +274,7 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 particleTypes = electrons
-runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux Ex
+runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir2d.py
 analysisOutputImage = langmuir2d_analysis.png
 tolerance = 1e-12
@@ -293,7 +293,7 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 particleTypes = electrons
-runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 warpx.do_dynamic_scheduling=0 diag1.fields_to_plot = Ex jx diag1.electrons.variables=w ux Ex
+runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 warpx.do_dynamic_scheduling=0 diag1.fields_to_plot = Ex jx diag1.electrons.variables=w ux
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir.py
 analysisOutputImage = langmuir_x_analysis.png
 tolerance = 5e-11
@@ -312,7 +312,7 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 particleTypes = electrons
-runtime_params = electrons.uy=0.01 electrons.ymax=0.e-6 warpx.do_dynamic_scheduling=0 diag1.fields_to_plot = Ey jy diag1.electrons.variables=w uy Ey
+runtime_params = electrons.uy=0.01 electrons.ymax=0.e-6 warpx.do_dynamic_scheduling=0 diag1.fields_to_plot = Ey jy diag1.electrons.variables=w uy
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir.py
 analysisOutputImage = langmuir_y_analysis.png
 tolerance = 5e-11
@@ -331,7 +331,7 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 particleTypes = electrons
-runtime_params = electrons.uz=0.01 electrons.zmax=0.e-6 warpx.do_dynamic_scheduling=0  diag1.fields_to_plot = Ez jz diag1.electrons.variables=w uz Ez
+runtime_params = electrons.uz=0.01 electrons.zmax=0.e-6 warpx.do_dynamic_scheduling=0  diag1.fields_to_plot = Ez jz diag1.electrons.variables=w uz
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir.py
 analysisOutputImage = langmuir_z_analysis.png
 tolerance = 5e-11
@@ -425,7 +425,7 @@ numthreads = 1
 compileTest = 0
 doVis = 0
 compareParticles = 0
-runtime_params = warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz Ex Ey Ez diag1.positrons.variables=w ux uy uz Ex Ey Ez
+runtime_params = warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz
 particleTypes = electrons positrons
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py
 analysisOutputImage = langmuir_multi_2d_analysis.png
@@ -444,7 +444,7 @@ numthreads = 1
 compileTest = 0
 doVis = 0
 compareParticles = 0
-runtime_params = psatd.fftw_plan_measure=0 diag1.electrons.variables=w ux uy uz Ex Ey Ez diag1.positrons.variables=w ux uy uz Ex Ey Ez diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell
+runtime_params = psatd.fftw_plan_measure=0 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
 particleTypes = electrons positrons
 analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py
 analysisOutputImage = langmuir_multi_2d_analysis.png
@@ -463,7 +463,7 @@ tolerance = 5e-11
 # compileTest = 0
 # doVis = 0
 # compareParticles = 0
-# runtime_params =  psatd.fftw_plan_measure=0 warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz Ex Ey Ez diag1.positrons.variables=w ux uy uz Ex Ey Ez diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell
+# runtime_params =  psatd.fftw_plan_measure=0 warpx.do_nodal=1 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
 # particleTypes = electrons positrons
 # analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_2d.py
 # analysisOutputImage = langmuir_multi_2d_analysis.png
@@ -481,7 +481,7 @@ tolerance = 5e-11
 # numthreads = 1
 # compileTest = 0
 # doVis = 0
-# runtime_params = diag1.electrons.variables=w ux uy uz Ex Ey Ez Bx By diag1.ions.variables=w ux uy uz Ex Ey Ez Bx By
+# runtime_params = diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz
 # compareParticles = 0
 # particleTypes = electrons ions
 # analysisRoutine = Examples/Tests/Langmuir/analysis_langmuir_multi_rz.py
@@ -623,7 +623,7 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 particleTypes = electrons
-outputFile = diags/plotfiles/plt00040
+outputFile = diags/diag200040
 tolerance = 5e-11
 
 [uniform_plasma_restart]
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index 9df25e1da..4789bc5fa 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -248,7 +248,7 @@ tolerance = 1.e-14
 [Langmuir_2d]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_rt
-runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux Ex
+runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux
 dim = 2
 addToCompileString =
 restartTest = 0
@@ -267,7 +267,7 @@ tolerance = 1.e-14
 [Langmuir_2d_single_precision]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_rt
-runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux Ex
+runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux
 dim = 2
 addToCompileString = PRECISION=FLOAT USE_SINGLE_PRECISION_PARTICLES=TRUE
 restartTest = 0
@@ -286,7 +286,7 @@ tolerance = 1.0e-4
 [Langmuir_2d_nompi]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_3d_rt
-runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux Ex
+runtime_params = electrons.ux=0.01 electrons.xmax=0.e-6 diag1.fields_to_plot=Ex jx diag1.electrons.variables=w ux
 dim = 2
 addToCompileString =
 restartTest = 0
@@ -457,7 +457,7 @@ tolerance = 5.e-7
 [Langmuir_multi_2d_nodal]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz Ex Ey Ez diag1.positrons.variables=w ux uy uz Ex Ey Ez
+runtime_params = warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz
 dim = 2
 addToCompileString =
 restartTest = 0
@@ -476,7 +476,7 @@ tolerance = 1.e-14
 [Langmuir_multi_2d_psatd]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = psatd.fftw_plan_measure=0 diag1.electrons.variables=w ux uy uz Ex Ey Ez diag1.positrons.variables=w ux uy uz Ex Ey Ez diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell
+runtime_params = psatd.fftw_plan_measure=0 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
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -495,7 +495,7 @@ tolerance = 1.e-14
 [Langmuir_multi_2d_psatd_momentum_conserving]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = algo.field_gathering=momentum-conserving psatd.fftw_plan_measure=0 diag1.electrons.variables=w ux uy uz Ex Ey Ez diag1.positrons.variables=w ux uy uz Ex Ey Ez diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell
+runtime_params = algo.field_gathering=momentum-conserving psatd.fftw_plan_measure=0 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
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -514,7 +514,7 @@ tolerance = 1.e-14
 [Langmuir_multi_2d_psatd_current_correction]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = amr.max_grid_size=128 psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.do_current_correction=1 diag1.electrons.variables=w ux uy uz Ex Ey Ez diag1.positrons.variables=w ux uy uz Ex Ey Ez diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE
+runtime_params = amr.max_grid_size=128 psatd.fftw_plan_measure=0 psatd.periodic_single_box_fft=1 psatd.do_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
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -532,7 +532,7 @@ analysisOutputImage = langmuir_multi_2d_analysis.png
 [Langmuir_multi_2d_psatd_nodal]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rt
-runtime_params = psatd.fftw_plan_measure=0 warpx.do_nodal=1 algo.current_deposition=direct diag1.electrons.variables=w ux uy uz Ex Ey Ez diag1.positrons.variables=w ux uy uz Ex Ey Ez diag1.fields_to_plot=Ex Ey Ez jx jy jz part_per_cell
+runtime_params = psatd.fftw_plan_measure=0 warpx.do_nodal=1 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
 dim = 2
 addToCompileString = USE_PSATD=TRUE
 restartTest = 0
@@ -551,7 +551,7 @@ tolerance = 1.e-14
 [Langmuir_multi_rz]
 buildDir = .
 inputFile = Examples/Tests/Langmuir/inputs_2d_multi_rz_rt
-runtime_params = diag1.electrons.variables=w ux uy uz Ex Ey Ez Bx By diag1.ions.variables=w ux uy uz Ex Ey Ez Bx By diag1.dump_rz_modes=0
+runtime_params = diag1.electrons.variables=w ux uy uz diag1.ions.variables=w ux uy uz diag1.dump_rz_modes=0
 dim = 2
 addToCompileString = USE_RZ=TRUE
 restartTest = 0
@@ -1105,7 +1105,7 @@ tolerance = 1.e-14
 [LaserAccelerationRZ]
 buildDir = .
 inputFile = Examples/Physics_applications/laser_acceleration/inputs_2d_rz
-runtime_params = diag1.electrons.variables=w ux uy uz Ex Ey Ez Bx By diag1.beam.variables=w ux uy uz Ex Ey Ez Bx By max_step=10
+runtime_params = diag1.electrons.variables=w ux uy uz diag1.beam.variables=w ux uy uz max_step=10
 dim = 2
 addToCompileString = USE_RZ=TRUE
 restartTest = 0
diff --git a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp
index fc685dc05..de548a848 100644
--- a/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp
+++ b/Source/Diagnostics/FlushFormats/FlushFormatPlotfile.cpp
@@ -275,14 +275,6 @@ FlushFormatPlotfile::WriteParticles(const std::string& dir,
         real_names.push_back("momentum_y");
         real_names.push_back("momentum_z");
 
-        real_names.push_back("Ex");
-        real_names.push_back("Ey");
-        real_names.push_back("Ez");
-
-        real_names.push_back("Bx");
-        real_names.push_back("By");
-        real_names.push_back("Bz");
-
 #ifdef WARPX_DIM_RZ
         real_names.push_back("theta");
 #endif
diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp
index 5615c09e8..fae5bb594 100644
--- a/Source/Diagnostics/FullDiagnostics.cpp
+++ b/Source/Diagnostics/FullDiagnostics.cpp
@@ -377,7 +377,4 @@ FullDiagnostics::PrepareFieldDataForOutput ()
     warpx.FillBoundaryAux(warpx.getngUpdateAux());
 #endif
     warpx.UpdateAuxilaryData();
-
-    warpx.FieldGather();
-
 }
diff --git a/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp b/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp
index 15bd01e86..b06817f91 100644
--- a/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp
+++ b/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp
@@ -10,7 +10,7 @@ ParticleDiag::ParticleDiag(std::string diag_name, std::string name, WarpXParticl
 {
     ParmParse pp(diag_name + "." + name);
     if (!pp.queryarr("variables", variables)){
-        variables = {"Ex", "Ey", "Ez", "Bx", "By", "Bz", "ux", "uy", "uz", "w"};
+        variables = {"ux", "uy", "uz", "w"};
     }
 
     //variable to set plot_flags size
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index 87d39ffc5..72621e52c 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -267,7 +267,8 @@ WarpX::OneStep_nosub (Real cur_time)
 
     // Loop over species. For each ionizable species, create particles in
     // product species.
-    mypc->doFieldIonization();
+    doFieldIonization();
+
     mypc->doCoulombCollisions();
 #ifdef WARPX_QED
     mypc->doQEDSchwinger();
@@ -298,8 +299,7 @@ WarpX::OneStep_nosub (Real cur_time)
 #endif
 
 #ifdef WARPX_QED
-    //Do QED processes
-    mypc->doQedEvents();
+    doQEDEvents();
 #endif
 
     // Synchronize J and rho
@@ -408,11 +408,10 @@ WarpX::OneStep_sub1 (Real curtime)
     // TODO: we could save some charge depositions
     // Loop over species. For each ionizable species, create particles in
     // product species.
-    mypc->doFieldIonization();
+    doFieldIonization();
 
 #ifdef WARPX_QED
-    //Do QED processes
-    mypc->doQedEvents();
+    doQEDEvents();
 #endif
 
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(finest_level == 1, "Must have exactly two levels");
@@ -550,6 +549,40 @@ WarpX::OneStep_sub1 (Real curtime)
 }
 
 void
+WarpX::doFieldIonization ()
+{
+    for (int lev = 0; lev <= finest_level; ++lev) {
+        doFieldIonization(lev);
+    }
+}
+
+void
+WarpX::doFieldIonization (int lev)
+{
+    mypc->doFieldIonization(lev,
+                            *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2],
+                            *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]);
+}
+
+#ifdef WARPX_QED
+void
+WarpX::doQEDEvents ()
+{
+    for (int lev = 0; lev <= finest_level; ++lev) {
+        doQEDEvents(lev);
+    }
+}
+
+void
+WarpX::doQEDEvents (int lev)
+{
+    mypc->doQedEvents(lev,
+                      *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2],
+                      *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]);
+}
+#endif
+
+void
 WarpX::PushParticlesandDepose (amrex::Real cur_time)
 {
     // Evolve particles to p^{n+1/2} and x^{n+1}
diff --git a/Source/Particles/CMakeLists.txt b/Source/Particles/CMakeLists.txt
index e4a665de5..0a8abe4eb 100644
--- a/Source/Particles/CMakeLists.txt
+++ b/Source/Particles/CMakeLists.txt
@@ -10,7 +10,7 @@ target_sources(WarpX
 add_subdirectory(Collision)
 #add_subdirectory(Deposition)
 add_subdirectory(ElementaryProcess)
-#add_subdirectory(Gather)
+add_subdirectory(Gather)
 add_subdirectory(ParticleCreation)
 #add_subdirectory(Pusher)
 add_subdirectory(Sorting)
diff --git a/Source/Particles/ElementaryProcess/CMakeLists.txt b/Source/Particles/ElementaryProcess/CMakeLists.txt
index 1a4550a94..bc3bff031 100644
--- a/Source/Particles/ElementaryProcess/CMakeLists.txt
+++ b/Source/Particles/ElementaryProcess/CMakeLists.txt
@@ -1,3 +1,13 @@
+target_sources(WarpX
+  PRIVATE
+    Ionization.cpp
+)
+
 if(WarpX_HAVE_QED)
+    target_sources(WarpX
+      PRIVATE
+        QEDPairGeneration.cpp
+        QEDPhotonEmission.cpp
+    )
     add_subdirectory(QEDInternals)
 endif()
diff --git a/Source/Particles/ElementaryProcess/Ionization.H b/Source/Particles/ElementaryProcess/Ionization.H
index 6cf30bd4d..6275c7af4 100644
--- a/Source/Particles/ElementaryProcess/Ionization.H
+++ b/Source/Particles/ElementaryProcess/Ionization.H
@@ -10,38 +10,96 @@
 
 #include "Utils/WarpXConst.H"
 #include "Particles/WarpXParticleContainer.H"
+#include "Particles/Gather/GetExternalFields.H"
+#include "Particles/Gather/FieldGather.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
 
 struct IonizationFilterFunc
 {
-    const amrex::Real* const AMREX_RESTRICT m_ionization_energies;
-    const amrex::Real* const AMREX_RESTRICT m_adk_prefactor;
-    const amrex::Real* const AMREX_RESTRICT m_adk_exp_prefactor;
-    const amrex::Real* const AMREX_RESTRICT m_adk_power;
+    const amrex::Real* AMREX_RESTRICT m_ionization_energies;
+    const amrex::Real* AMREX_RESTRICT m_adk_prefactor;
+    const amrex::Real* AMREX_RESTRICT m_adk_exp_prefactor;
+    const amrex::Real* AMREX_RESTRICT m_adk_power;
 
     int comp;
     int m_atomic_number;
 
+    GetParticlePosition m_get_position;
+    GetExternalEField m_get_externalE;
+    GetExternalBField m_get_externalB;
+
+    amrex::Array4<const amrex::Real> m_ex_arr;
+    amrex::Array4<const amrex::Real> m_ey_arr;
+    amrex::Array4<const amrex::Real> m_ez_arr;
+    amrex::Array4<const amrex::Real> m_bx_arr;
+    amrex::Array4<const amrex::Real> m_by_arr;
+    amrex::Array4<const amrex::Real> m_bz_arr;
+
+    amrex::IndexType m_ex_type;
+    amrex::IndexType m_ey_type;
+    amrex::IndexType m_ez_type;
+    amrex::IndexType m_bx_type;
+    amrex::IndexType m_by_type;
+    amrex::IndexType m_bz_type;
+
+    amrex::GpuArray<amrex::Real, 3> m_dx_arr;
+    amrex::GpuArray<amrex::Real, 3> m_xyzmin_arr;
+
+    int m_l_lower_order_in_v;
+    int m_nox;
+    int m_n_rz_azimuthal_modes;
+
+    amrex::Dim3 m_lo;
+
+    IonizationFilterFunc (const WarpXParIter& a_pti, int lev, int ngE,
+                          amrex::FArrayBox const& exfab,
+                          amrex::FArrayBox const& eyfab,
+                          amrex::FArrayBox const& ezfab,
+                          amrex::FArrayBox const& bxfab,
+                          amrex::FArrayBox const& byfab,
+                          amrex::FArrayBox const& bzfab,
+                          amrex::Array<amrex::Real,3> v_galilean,
+                          const amrex::Real* const AMREX_RESTRICT a_ionization_energies,
+                          const amrex::Real* const AMREX_RESTRICT a_adk_prefactor,
+                          const amrex::Real* const AMREX_RESTRICT a_adk_exp_prefactor,
+                          const amrex::Real* const AMREX_RESTRICT a_adk_power,
+                          int a_comp,
+                          int a_atomic_number,
+                          int a_offset = 0) noexcept;
+
     template <typename PData>
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     bool operator() (const PData& ptd, int i) const noexcept
     {
+        using namespace amrex::literals;
+
         const int ion_lev = ptd.m_runtime_idata[comp][i];
         if (ion_lev < m_atomic_number)
         {
             constexpr amrex::Real c = PhysConst::c;
             constexpr amrex::Real c2_inv = 1./c/c;
 
+            // gather E and B
+            amrex::ParticleReal xp, yp, zp;
+            m_get_position(i, xp, yp, zp);
+
+            amrex::ParticleReal ex = 0._rt, ey = 0._rt, ez = 0._rt;
+            m_get_externalE(i, ex, ey, ez);
+
+            amrex::ParticleReal bx = 0._rt, by = 0._rt, bz = 0._rt;
+            m_get_externalB(i, bx, by, bz);
+
+            doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz,
+                           m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr,
+                           m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type,
+                           m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes,
+                           m_nox, m_l_lower_order_in_v);
+
             // Compute electric field amplitude in the particle's frame of
             // reference (particularly important when in boosted frame).
             amrex::ParticleReal ux = ptd.m_rdata[PIdx::ux][i];
             amrex::ParticleReal uy = ptd.m_rdata[PIdx::uy][i];
             amrex::ParticleReal uz = ptd.m_rdata[PIdx::uz][i];
-            amrex::ParticleReal ex = ptd.m_rdata[PIdx::Ex][i];
-            amrex::ParticleReal ey = ptd.m_rdata[PIdx::Ey][i];
-            amrex::ParticleReal ez = ptd.m_rdata[PIdx::Ez][i];
-            amrex::ParticleReal bx = ptd.m_rdata[PIdx::Bx][i];
-            amrex::ParticleReal by = ptd.m_rdata[PIdx::By][i];
-            amrex::ParticleReal bz = ptd.m_rdata[PIdx::Bz][i];
 
             amrex::Real ga = std::sqrt(1. + (ux*ux + uy*uy + uz*uz) * c2_inv);
             amrex::Real E = std::sqrt(
diff --git a/Source/Particles/ElementaryProcess/Ionization.cpp b/Source/Particles/ElementaryProcess/Ionization.cpp
new file mode 100644
index 000000000..4c30987f4
--- /dev/null
+++ b/Source/Particles/ElementaryProcess/Ionization.cpp
@@ -0,0 +1,72 @@
+/* Copyright 2019-2020 Andrew Myers, Axel Huebl,
+ * Maxence Thevenet
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "WarpX.H"
+#include "Particles/ElementaryProcess/Ionization.H"
+
+IonizationFilterFunc::IonizationFilterFunc (const WarpXParIter& a_pti, int lev, int ngE,
+                                            amrex::FArrayBox const& exfab,
+                                            amrex::FArrayBox const& eyfab,
+                                            amrex::FArrayBox const& ezfab,
+                                            amrex::FArrayBox const& bxfab,
+                                            amrex::FArrayBox const& byfab,
+                                            amrex::FArrayBox const& bzfab,
+                                            amrex::Array<amrex::Real,3> v_galilean,
+                                            const amrex::Real* const AMREX_RESTRICT a_ionization_energies,
+                                            const amrex::Real* const AMREX_RESTRICT a_adk_prefactor,
+                                            const amrex::Real* const AMREX_RESTRICT a_adk_exp_prefactor,
+                                            const amrex::Real* const AMREX_RESTRICT a_adk_power,
+                                            int a_comp,
+                                            int a_atomic_number,
+                                            int a_offset) noexcept
+{
+    m_ionization_energies = a_ionization_energies;
+    m_adk_prefactor = a_adk_prefactor;
+    m_adk_exp_prefactor = a_adk_exp_prefactor;
+    m_adk_power = a_adk_power;
+    comp = a_comp;
+    m_atomic_number = a_atomic_number;
+
+    m_get_position  = GetParticlePosition(a_pti, a_offset);
+    m_get_externalE = GetExternalEField  (a_pti, a_offset);
+    m_get_externalB = GetExternalBField  (a_pti, a_offset);
+
+    m_ex_arr = exfab.array();
+    m_ey_arr = eyfab.array();
+    m_ez_arr = ezfab.array();
+    m_bx_arr = bxfab.array();
+    m_by_arr = byfab.array();
+    m_bz_arr = bzfab.array();
+
+    m_ex_type = exfab.box().ixType();
+    m_ey_type = eyfab.box().ixType();
+    m_ez_type = ezfab.box().ixType();
+    m_bx_type = bxfab.box().ixType();
+    m_by_type = byfab.box().ixType();
+    m_bz_type = bzfab.box().ixType();
+
+    amrex::Box box = a_pti.tilebox();
+    box.grow(ngE);
+
+    const std::array<amrex::Real,3>& dx = WarpX::CellSize(std::max(lev, 0));
+    m_dx_arr = {dx[0], dx[1], dx[2]};
+
+    // Lower corner of tile box physical domain (take into account Galilean shift)
+    amrex::Real cur_time = WarpX::GetInstance().gett_new(lev);
+    const auto& time_of_last_gal_shift = WarpX::GetInstance().time_of_last_gal_shift;
+    amrex::Real time_shift = (cur_time - time_of_last_gal_shift);
+    amrex::Array<amrex::Real,3> galilean_shift = { v_galilean[0]*time_shift, v_galilean[1]*time_shift, v_galilean[2]*time_shift };
+    const std::array<amrex::Real, 3>& xyzmin = WarpX::LowerCorner(box, galilean_shift, lev);
+    m_xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+    m_l_lower_order_in_v = WarpX::l_lower_order_in_v;
+    m_nox = WarpX::nox;
+    m_n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+    m_lo = amrex::lbound(box);
+}
diff --git a/Source/Particles/ElementaryProcess/Make.package b/Source/Particles/ElementaryProcess/Make.package
index 588d017e6..cc0610731 100644
--- a/Source/Particles/ElementaryProcess/Make.package
+++ b/Source/Particles/ElementaryProcess/Make.package
@@ -1,4 +1,8 @@
+CEXE_sources += Ionization.cpp
+
 ifeq ($(QED),TRUE)
+     CEXE_sources += QEDPairGeneration.cpp
+     CEXE_sources += QEDPhotonEmission.cpp
     include $(WARPX_HOME)/Source/Particles/ElementaryProcess/QEDInternals/Make.package
     VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Particles/ElementaryProcess/QEDInternals/
 endif
diff --git a/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H b/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H
index bed6337c3..7056ac884 100644
--- a/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H
+++ b/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H
@@ -84,6 +84,9 @@ public:
 class BreitWheelerEvolveOpticalDepth
 {
 public:
+
+    BreitWheelerEvolveOpticalDepth () = default;
+
     /**
      * Constructor acquires a reference to control parameters and
      * lookup tables data.
@@ -138,10 +141,10 @@ public:
     }
 
 private:
-    //laser wavelenght is not used with SI units
-    const amrex::Real m_dummy_lambda{1.0};
+    //laser wavelength is not used with SI units
+    amrex::Real m_dummy_lambda{1.0};
 
-    const PicsarBreitWheelerCtrl m_ctrl;
+    PicsarBreitWheelerCtrl m_ctrl;
 
     //lookup table data
     size_t m_TTfunc_size;
diff --git a/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H b/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H
index 6d46c9019..fe7d560cd 100644
--- a/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H
+++ b/Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H
@@ -84,6 +84,9 @@ public:
 class QuantumSynchrotronEvolveOpticalDepth
 {
 public:
+
+    QuantumSynchrotronEvolveOpticalDepth () = default;
+
     /**
      * Constructor acquires pointers to control parameters and
      * lookup tables data.
@@ -139,10 +142,10 @@ public:
     }
 
 private:
-    //laser wavelenght is not used with SI units
-    const amrex::Real m_dummy_lambda{1.0};
+    //laser wavelength is not used with SI units
+    amrex::Real m_dummy_lambda{1.0};
 
-    const PicsarQuantumSynchrotronCtrl m_ctrl;
+    PicsarQuantumSynchrotronCtrl m_ctrl;
 
     //lookup table data
     size_t m_KKfunc_size;
diff --git a/Source/Particles/ElementaryProcess/QEDPairGeneration.H b/Source/Particles/ElementaryProcess/QEDPairGeneration.H
index 22f37a351..838afcaac 100644
--- a/Source/Particles/ElementaryProcess/QEDPairGeneration.H
+++ b/Source/Particles/ElementaryProcess/QEDPairGeneration.H
@@ -10,7 +10,9 @@
 
 #include "Utils/WarpXConst.H"
 #include "Particles/WarpXParticleContainer.H"
-
+#include "Particles/Gather/GetExternalFields.H"
+#include "Particles/Gather/FieldGather.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
 #include "QEDInternals/BreitWheelerEngineWrapper.H"
 
 /** @file
@@ -32,9 +34,9 @@ public:
     *
     * @param[in] opt_depth_runtime_comp index of the optical depth component
     */
-    PairGenerationFilterFunc(int const opt_depth_runtime_comp):
-        m_opt_depth_runtime_comp{opt_depth_runtime_comp}
-    {}
+    PairGenerationFilterFunc(int const opt_depth_runtime_comp)
+        : m_opt_depth_runtime_comp(opt_depth_runtime_comp)
+        {}
 
     /**
     * \brief Functor call. This method determines if a given (photon) particle
@@ -75,9 +77,16 @@ public:
     *
     * @param[in] generate_functor functor to be called to determine the properties of the generated pairs
     */
-    PairGenerationTransformFunc(BreitWheelerGeneratePairs const generate_functor):
-        m_generate_functor{generate_functor}
-    {}
+    PairGenerationTransformFunc(BreitWheelerGeneratePairs const generate_functor,
+                                const WarpXParIter& a_pti, int lev, int ngE,
+                                amrex::FArrayBox const& exfab,
+                                amrex::FArrayBox const& eyfab,
+                                amrex::FArrayBox const& ezfab,
+                                amrex::FArrayBox const& bxfab,
+                                amrex::FArrayBox const& byfab,
+                                amrex::FArrayBox const& bzfab,
+                                amrex::Array<amrex::Real,3> v_galilean,
+                                int a_offset = 0);
 
     /**
     * \brief Functor call. It determines the properties of the generated pair
@@ -104,12 +113,22 @@ public:
         const ParticleReal ux = src.m_rdata[PIdx::ux][i_src];
         const ParticleReal uy = src.m_rdata[PIdx::uy][i_src];
         const ParticleReal uz = src.m_rdata[PIdx::uz][i_src];
-        const ParticleReal ex = src.m_rdata[PIdx::Ex][i_src];
-        const ParticleReal ey = src.m_rdata[PIdx::Ey][i_src];
-        const ParticleReal ez = src.m_rdata[PIdx::Ez][i_src];
-        const ParticleReal bx = src.m_rdata[PIdx::Bx][i_src];
-        const ParticleReal by = src.m_rdata[PIdx::By][i_src];
-        const ParticleReal bz = src.m_rdata[PIdx::Bz][i_src];
+
+        // gather E and B
+        amrex::ParticleReal xp, yp, zp;
+        m_get_position(i_src, xp, yp, zp);
+
+        amrex::ParticleReal ex = 0._rt, ey = 0._rt, ez = 0._rt;
+        m_get_externalE(i_src, ex, ey, ez);
+
+        amrex::ParticleReal bx = 0._rt, by = 0._rt, bz = 0._rt;
+        m_get_externalB(i_src, bx, by, bz);
+
+        doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz,
+                       m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr,
+                       m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type,
+                       m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes,
+                       m_nox, m_l_lower_order_in_v);
 
         const auto px = ux*me;
         const auto py = uy*me;
@@ -153,7 +172,34 @@ public:
 private:
 
     const BreitWheelerGeneratePairs
-        m_generate_functor; /*!< A copy of the functor to generate pairs. It contains only pointers to the lookup tables.*/
+    m_generate_functor; /*!< A copy of the functor to generate pairs. It contains only pointers to the lookup tables.*/
+
+    GetParticlePosition m_get_position;
+    GetExternalEField m_get_externalE;
+    GetExternalBField m_get_externalB;
+
+    amrex::Array4<const amrex::Real> m_ex_arr;
+    amrex::Array4<const amrex::Real> m_ey_arr;
+    amrex::Array4<const amrex::Real> m_ez_arr;
+    amrex::Array4<const amrex::Real> m_bx_arr;
+    amrex::Array4<const amrex::Real> m_by_arr;
+    amrex::Array4<const amrex::Real> m_bz_arr;
+
+    amrex::IndexType m_ex_type;
+    amrex::IndexType m_ey_type;
+    amrex::IndexType m_ez_type;
+    amrex::IndexType m_bx_type;
+    amrex::IndexType m_by_type;
+    amrex::IndexType m_bz_type;
+
+    amrex::GpuArray<amrex::Real, 3> m_dx_arr;
+    amrex::GpuArray<amrex::Real, 3> m_xyzmin_arr;
+
+    int m_l_lower_order_in_v;
+    int m_nox;
+    int m_n_rz_azimuthal_modes;
+
+    amrex::Dim3 m_lo;
 };
 
 #endif //QED_PAIR_GENERATION_H_
diff --git a/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp b/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp
new file mode 100644
index 000000000..cdb88d281
--- /dev/null
+++ b/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp
@@ -0,0 +1,62 @@
+/* Copyright 2019-2020 Andrew Myers, Axel Huebl,
+ * Maxence Thevenet
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "WarpX.H"
+#include "Particles/ElementaryProcess/QEDPairGeneration.H"
+
+PairGenerationTransformFunc::
+PairGenerationTransformFunc (BreitWheelerGeneratePairs const generate_functor,
+                             const WarpXParIter& a_pti, int lev, int ngE,
+                             amrex::FArrayBox const& exfab,
+                             amrex::FArrayBox const& eyfab,
+                             amrex::FArrayBox const& ezfab,
+                             amrex::FArrayBox const& bxfab,
+                             amrex::FArrayBox const& byfab,
+                             amrex::FArrayBox const& bzfab,
+                             amrex::Array<amrex::Real,3> v_galilean,
+                             int a_offset)
+: m_generate_functor(generate_functor)
+{
+    m_get_position  = GetParticlePosition(a_pti, a_offset);
+    m_get_externalE = GetExternalEField  (a_pti, a_offset);
+    m_get_externalB = GetExternalBField  (a_pti, a_offset);
+
+    m_ex_arr = exfab.array();
+    m_ey_arr = eyfab.array();
+    m_ez_arr = ezfab.array();
+    m_bx_arr = bxfab.array();
+    m_by_arr = byfab.array();
+    m_bz_arr = bzfab.array();
+
+    m_ex_type = exfab.box().ixType();
+    m_ey_type = eyfab.box().ixType();
+    m_ez_type = ezfab.box().ixType();
+    m_bx_type = bxfab.box().ixType();
+    m_by_type = byfab.box().ixType();
+    m_bz_type = bzfab.box().ixType();
+
+    amrex::Box box = a_pti.tilebox();
+    box.grow(ngE);
+
+    const std::array<amrex::Real,3>& dx = WarpX::CellSize(std::max(lev, 0));
+    m_dx_arr = {dx[0], dx[1], dx[2]};
+
+    // Lower corner of tile box physical domain (take into account Galilean shift)
+    amrex::Real cur_time = WarpX::GetInstance().gett_new(lev);
+    const auto& time_of_last_gal_shift = WarpX::GetInstance().time_of_last_gal_shift;
+    amrex::Real time_shift = (cur_time - time_of_last_gal_shift);
+    amrex::Array<amrex::Real,3> galilean_shift = { v_galilean[0]*time_shift, v_galilean[1]*time_shift, v_galilean[2]*time_shift };
+    const std::array<amrex::Real, 3>& xyzmin = WarpX::LowerCorner(box, galilean_shift, lev);
+    m_xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+    m_l_lower_order_in_v = WarpX::l_lower_order_in_v;
+    m_nox = WarpX::nox;
+    m_n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+    m_lo = amrex::lbound(box);
+}
diff --git a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H
index c7c801eeb..4dc573f60 100644
--- a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H
+++ b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H
@@ -10,7 +10,9 @@
 
 #include "Utils/WarpXConst.H"
 #include "Particles/WarpXParticleContainer.H"
-
+#include "Particles/Gather/GetExternalFields.H"
+#include "Particles/Gather/FieldGather.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
 #include "QEDInternals/QuantumSyncEngineWrapper.H"
 
 /** @file
@@ -32,9 +34,9 @@ public:
     *
     * @param[in] opt_depth_runtime_comp Index of the optical depth component
     */
-    PhotonEmissionFilterFunc(int const opt_depth_runtime_comp):
-        m_opt_depth_runtime_comp{opt_depth_runtime_comp}
-    {}
+    PhotonEmissionFilterFunc(int const opt_depth_runtime_comp)
+        : m_opt_depth_runtime_comp(opt_depth_runtime_comp)
+        {}
 
     /**
     * \brief Functor call. This method determines if a given (electron or positron)
@@ -56,7 +58,7 @@ public:
     }
 
 private:
-    int m_opt_depth_runtime_comp /*!< Index of the optical depth component of the source species*/;
+    int m_opt_depth_runtime_comp; /*!< Index of the optical depth component of the source species*/
 };
 
 /**
@@ -81,15 +83,19 @@ public:
     * @param[in] opt_depth_runtime_comp index of the optical depth component of the source species
     * @param[in] emission_functor functor to generate photons and update momentum of the source particles
     */
-    PhotonEmissionTransformFunc(
+    PhotonEmissionTransformFunc (
         QuantumSynchrotronGetOpticalDepth opt_depth_functor,
         int const opt_depth_runtime_comp,
-        QuantumSynchrotronGeneratePhotonAndUpdateMomentum const emission_functor
-        ):
-        m_opt_depth_functor{opt_depth_functor},
-        m_opt_depth_runtime_comp{opt_depth_runtime_comp},
-        m_emission_functor{emission_functor}
-        {}
+        QuantumSynchrotronGeneratePhotonAndUpdateMomentum const emission_functor,
+        const WarpXParIter& a_pti, int lev, int ngE,
+        amrex::FArrayBox const& exfab,
+        amrex::FArrayBox const& eyfab,
+        amrex::FArrayBox const& ezfab,
+        amrex::FArrayBox const& bxfab,
+        amrex::FArrayBox const& byfab,
+        amrex::FArrayBox const& bzfab,
+        amrex::Array<amrex::Real,3> v_galilean,
+        int a_offset = 0);
 
     /**
     * \brief Functor call. It determines the properties of the generated photon
@@ -113,12 +119,22 @@ public:
         const ParticleReal ux = src.m_rdata[PIdx::ux][i_src];
         const ParticleReal uy = src.m_rdata[PIdx::uy][i_src];
         const ParticleReal uz = src.m_rdata[PIdx::uz][i_src];
-        const ParticleReal ex = src.m_rdata[PIdx::Ex][i_src];
-        const ParticleReal ey = src.m_rdata[PIdx::Ey][i_src];
-        const ParticleReal ez = src.m_rdata[PIdx::Ez][i_src];
-        const ParticleReal bx = src.m_rdata[PIdx::Bx][i_src];
-        const ParticleReal by = src.m_rdata[PIdx::By][i_src];
-        const ParticleReal bz = src.m_rdata[PIdx::Bz][i_src];
+
+        // gather E and B
+        amrex::ParticleReal xp, yp, zp;
+        m_get_position(i_src, xp, yp, zp);
+
+        amrex::ParticleReal ex = 0._rt, ey = 0._rt, ez = 0._rt;
+        m_get_externalE(i_src, ex, ey, ez);
+
+        amrex::ParticleReal bx = 0._rt, by = 0._rt, bz = 0._rt;
+        m_get_externalB(i_src, bx, by, bz);
+
+        doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz,
+                       m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr,
+                       m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type,
+                       m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes,
+                       m_nox, m_l_lower_order_in_v);
 
         // Particle momentum is stored as gamma * velocity.
         // Convert to m * gamma * velocity before applying the emission functor.
@@ -160,6 +176,33 @@ private:
     const QuantumSynchrotronGeneratePhotonAndUpdateMomentum
         m_emission_functor;  /*!< A copy of the functor to generate photons. It contains only pointers to the lookup tables.*/
     const int m_opt_depth_runtime_comp = 0;  /*!< Index of the optical depth component of source species*/
+
+    GetParticlePosition m_get_position;
+    GetExternalEField m_get_externalE;
+    GetExternalBField m_get_externalB;
+
+    amrex::Array4<const amrex::Real> m_ex_arr;
+    amrex::Array4<const amrex::Real> m_ey_arr;
+    amrex::Array4<const amrex::Real> m_ez_arr;
+    amrex::Array4<const amrex::Real> m_bx_arr;
+    amrex::Array4<const amrex::Real> m_by_arr;
+    amrex::Array4<const amrex::Real> m_bz_arr;
+
+    amrex::IndexType m_ex_type;
+    amrex::IndexType m_ey_type;
+    amrex::IndexType m_ez_type;
+    amrex::IndexType m_bx_type;
+    amrex::IndexType m_by_type;
+    amrex::IndexType m_bz_type;
+
+    amrex::GpuArray<amrex::Real, 3> m_dx_arr;
+    amrex::GpuArray<amrex::Real, 3> m_xyzmin_arr;
+
+    int m_l_lower_order_in_v;
+    int m_nox;
+    int m_n_rz_azimuthal_modes;
+
+    amrex::Dim3 m_lo;
 };
 
 
diff --git a/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp b/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp
new file mode 100644
index 000000000..6134ff542
--- /dev/null
+++ b/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp
@@ -0,0 +1,66 @@
+/* Copyright 2019-2020 Andrew Myers, Axel Huebl,
+ * Maxence Thevenet
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "WarpX.H"
+#include "Particles/ElementaryProcess/QEDPhotonEmission.H"
+
+PhotonEmissionTransformFunc::
+PhotonEmissionTransformFunc (QuantumSynchrotronGetOpticalDepth opt_depth_functor,
+                             int const opt_depth_runtime_comp,
+                             QuantumSynchrotronGeneratePhotonAndUpdateMomentum const emission_functor,
+                             const WarpXParIter& a_pti, int lev, int ngE,
+                             amrex::FArrayBox const& exfab,
+                             amrex::FArrayBox const& eyfab,
+                             amrex::FArrayBox const& ezfab,
+                             amrex::FArrayBox const& bxfab,
+                             amrex::FArrayBox const& byfab,
+                             amrex::FArrayBox const& bzfab,
+                             amrex::Array<amrex::Real,3> v_galilean,
+                             int a_offset)
+:m_opt_depth_functor{opt_depth_functor},
+ m_opt_depth_runtime_comp{opt_depth_runtime_comp},
+ m_emission_functor{emission_functor}
+{
+    m_get_position  = GetParticlePosition(a_pti, a_offset);
+    m_get_externalE = GetExternalEField  (a_pti, a_offset);
+    m_get_externalB = GetExternalBField  (a_pti, a_offset);
+
+    m_ex_arr = exfab.array();
+    m_ey_arr = eyfab.array();
+    m_ez_arr = ezfab.array();
+    m_bx_arr = bxfab.array();
+    m_by_arr = byfab.array();
+    m_bz_arr = bzfab.array();
+
+    m_ex_type = exfab.box().ixType();
+    m_ey_type = eyfab.box().ixType();
+    m_ez_type = ezfab.box().ixType();
+    m_bx_type = bxfab.box().ixType();
+    m_by_type = byfab.box().ixType();
+    m_bz_type = bzfab.box().ixType();
+
+    amrex::Box box = a_pti.tilebox();
+    box.grow(ngE);
+
+    const std::array<amrex::Real,3>& dx = WarpX::CellSize(std::max(lev, 0));
+    m_dx_arr = {dx[0], dx[1], dx[2]};
+
+    // Lower corner of tile box physical domain (take into account Galilean shift)
+    amrex::Real cur_time = WarpX::GetInstance().gett_new(lev);
+    const auto& time_of_last_gal_shift = WarpX::GetInstance().time_of_last_gal_shift;
+    amrex::Real time_shift = (cur_time - time_of_last_gal_shift);
+    amrex::Array<amrex::Real,3> galilean_shift = { v_galilean[0]*time_shift, v_galilean[1]*time_shift, v_galilean[2]*time_shift };
+    const std::array<amrex::Real, 3>& xyzmin = WarpX::LowerCorner(box, galilean_shift, lev);
+    m_xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+    m_l_lower_order_in_v = WarpX::l_lower_order_in_v;
+    m_nox = WarpX::nox;
+    m_n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+    m_lo = amrex::lbound(box);
+}
diff --git a/Source/Particles/Gather/CMakeLists.txt b/Source/Particles/Gather/CMakeLists.txt
new file mode 100644
index 000000000..8ede378f2
--- /dev/null
+++ b/Source/Particles/Gather/CMakeLists.txt
@@ -0,0 +1,4 @@
+target_sources(WarpX
+  PRIVATE
+    GetExternalFields.cpp
+)
diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H
index fa92a7ae9..551c96ac1 100644
--- a/Source/Particles/Gather/FieldGather.H
+++ b/Source/Particles/Gather/FieldGather.H
@@ -16,7 +16,7 @@
 /**
  * \brief Field gather for a single particle
  *
- * /param xp, yp, zp                : Particle position coordinates
+ * \param xp, yp, zp                : Particle position coordinates
  * \param Exp, Eyp, Ezp             : Electric field on particles.
  * \param Bxp, Byp, Bzp             : Magnetic field on particles.
  * \param ex_arr ey_arr ez_arr      : Array4 of the electric field, either full array or tile.
@@ -430,4 +430,87 @@ void doGatherShapeN(const GetParticlePosition& getPosition,
         );
 }
 
+/**
+ * \brief Field gather for a single particle
+ *
+ * /param xp, yp, zp                : Particle position coordinates
+ * \param Exp, Eyp, Ezp             : Electric field on particles.
+ * \param Bxp, Byp, Bzp             : Magnetic field on particles.
+ * \param ex_arr ey_arr ez_arr      : Array4 of the electric field, either full array or tile.
+ * \param bx_arr by_arr bz_arr      : Array4 of the magnetic field, either full array or tile.
+ * \param ex_type, ey_type, ez_type : IndexType of the electric field
+ * \param bx_type, by_type, bz_type : IndexType of the magnetic field
+ * \param dx                        : 3D cell spacing
+ * \param xyzmin                    : Physical lower bounds of domain in x, y, z.
+ * \param lo                        : Index lower bounds of domain.
+ * \param n_rz_azimuthal_modes      : Number of azimuthal modes when using RZ geometry
+ * \param nox                       : order of the particle shape function
+ * \param l_lower_order_in_v        : whether to use lower order in v
+ */
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void doGatherShapeN (const amrex::ParticleReal xp,
+                     const amrex::ParticleReal yp,
+                     const amrex::ParticleReal zp,
+                     amrex::ParticleReal& Exp,
+                     amrex::ParticleReal& Eyp,
+                     amrex::ParticleReal& Ezp,
+                     amrex::ParticleReal& Bxp,
+                     amrex::ParticleReal& Byp,
+                     amrex::ParticleReal& Bzp,
+                     amrex::Array4<amrex::Real const> const& ex_arr,
+                     amrex::Array4<amrex::Real const> const& ey_arr,
+                     amrex::Array4<amrex::Real const> const& ez_arr,
+                     amrex::Array4<amrex::Real const> const& bx_arr,
+                     amrex::Array4<amrex::Real const> const& by_arr,
+                     amrex::Array4<amrex::Real const> const& bz_arr,
+                     const amrex::IndexType ex_type,
+                     const amrex::IndexType ey_type,
+                     const amrex::IndexType ez_type,
+                     const amrex::IndexType bx_type,
+                     const amrex::IndexType by_type,
+                     const amrex::IndexType bz_type,
+                     const amrex::GpuArray<amrex::Real, 3>& dx_arr,
+                     const amrex::GpuArray<amrex::Real, 3>& xyzmin_arr,
+                     const amrex::Dim3& lo,
+                     const long n_rz_azimuthal_modes,
+                     const int nox,
+                     const int l_lower_order_in_v)
+{
+    if (l_lower_order_in_v) {
+        if (nox == 1) {
+            doGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                                ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                                ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                                dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
+        } else if (nox == 2) {
+            doGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                                ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                                ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                                dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
+        } else if (nox == 3) {
+            doGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                                ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                                ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                                dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
+        }
+    } else {
+        if (nox == 1) {
+            doGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                                ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                                ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                                dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
+        } else if (nox == 2) {
+            doGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                                ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                                ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                                dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
+        } else if (nox == 3) {
+            doGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                                ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                                ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                                dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes);
+        }
+    }
+}
+
 #endif // FIELDGATHER_H_
diff --git a/Source/Particles/Gather/GetExternalFields.H b/Source/Particles/Gather/GetExternalFields.H
index 8e08fc101..f762febd5 100644
--- a/Source/Particles/Gather/GetExternalFields.H
+++ b/Source/Particles/Gather/GetExternalFields.H
@@ -2,6 +2,7 @@
 #define WARPX_PARTICLES_GATHER_GETEXTERNALFIELDS_H_
 
 #include "Particles/WarpXParticleContainer.H"
+#include "Particles/Pusher/GetAndSetPosition.H"
 
 #include <AMReX_REAL.H>
 
@@ -60,27 +61,9 @@ struct GetExternalField
 */
 struct GetExternalEField : GetExternalField
 {
-    GetExternalEField (const WarpXParIter& a_pti, int a_offset = 0) noexcept
-    {
-        auto& warpx = WarpX::GetInstance();
-        auto& mypc = warpx.GetPartContainer();
-        if (mypc.m_E_ext_particle_s=="constant" || mypc.m_E_ext_particle_s=="default")
-        {
-            m_type = Constant;
-            m_field_value[0] = mypc.m_E_external_particle[0];
-            m_field_value[1] = mypc.m_E_external_particle[1];
-            m_field_value[2] = mypc.m_E_external_particle[2];
-        }
-        else if (mypc.m_E_ext_particle_s=="parse_e_ext_particle_function")
-        {
-            m_type = Parser;
-            m_time = warpx.gett_new(a_pti.GetLevel());
-            m_get_position = GetParticlePosition(a_pti, a_offset);
-            m_xfield_partparser = mypc.m_Ex_particle_parser.get();
-            m_yfield_partparser = mypc.m_Ey_particle_parser.get();
-            m_zfield_partparser = mypc.m_Ez_particle_parser.get();
-        }
-    }
+    GetExternalEField () = default;
+
+    GetExternalEField (const WarpXParIter& a_pti, int a_offset = 0) noexcept;
 };
 
 /** \brief Functor that can be used to assign the external
@@ -88,27 +71,9 @@ struct GetExternalEField : GetExternalField
 */
 struct GetExternalBField : GetExternalField
 {
-    GetExternalBField (const WarpXParIter& a_pti, int a_offset = 0) noexcept
-    {
-        auto& warpx = WarpX::GetInstance();
-        auto& mypc = warpx.GetPartContainer();
-        if (mypc.m_B_ext_particle_s=="constant" || mypc.m_B_ext_particle_s=="default")
-        {
-            m_type = Constant;
-            m_field_value[0] = mypc.m_B_external_particle[0];
-            m_field_value[1] = mypc.m_B_external_particle[1];
-            m_field_value[2] = mypc.m_B_external_particle[2];
-        }
-        else if (mypc.m_B_ext_particle_s=="parse_e_ext_particle_function")
-        {
-            m_type = Parser;
-            m_time = warpx.gett_new(a_pti.GetLevel());
-            m_get_position = GetParticlePosition(a_pti, a_offset);
-            m_xfield_partparser = mypc.m_Bx_particle_parser.get();
-            m_yfield_partparser = mypc.m_By_particle_parser.get();
-            m_zfield_partparser = mypc.m_Bz_particle_parser.get();
-        }
-    }
+    GetExternalBField () = default;
+
+    GetExternalBField (const WarpXParIter& a_pti, int a_offset = 0) noexcept;
 };
 
 #endif
diff --git a/Source/Particles/Gather/GetExternalFields.cpp b/Source/Particles/Gather/GetExternalFields.cpp
new file mode 100644
index 000000000..6e4eb6d1b
--- /dev/null
+++ b/Source/Particles/Gather/GetExternalFields.cpp
@@ -0,0 +1,46 @@
+#include "WarpX.H"
+#include "Particles/Gather/GetExternalFields.H"
+
+GetExternalEField::GetExternalEField (const WarpXParIter& a_pti, int a_offset) noexcept
+{
+    auto& warpx = WarpX::GetInstance();
+    auto& mypc = warpx.GetPartContainer();
+    if (mypc.m_E_ext_particle_s=="constant" || mypc.m_E_ext_particle_s=="default")
+    {
+        m_type = Constant;
+        m_field_value[0] = mypc.m_E_external_particle[0];
+        m_field_value[1] = mypc.m_E_external_particle[1];
+        m_field_value[2] = mypc.m_E_external_particle[2];
+    }
+    else if (mypc.m_E_ext_particle_s=="parse_e_ext_particle_function")
+    {
+        m_type = Parser;
+        m_time = warpx.gett_new(a_pti.GetLevel());
+        m_get_position = GetParticlePosition(a_pti, a_offset);
+        m_xfield_partparser = mypc.m_Ex_particle_parser.get();
+        m_yfield_partparser = mypc.m_Ey_particle_parser.get();
+        m_zfield_partparser = mypc.m_Ez_particle_parser.get();
+    }
+}
+
+GetExternalBField::GetExternalBField (const WarpXParIter& a_pti, int a_offset) noexcept
+{
+    auto& warpx = WarpX::GetInstance();
+    auto& mypc = warpx.GetPartContainer();
+    if (mypc.m_B_ext_particle_s=="constant" || mypc.m_B_ext_particle_s=="default")
+    {
+        m_type = Constant;
+        m_field_value[0] = mypc.m_B_external_particle[0];
+        m_field_value[1] = mypc.m_B_external_particle[1];
+        m_field_value[2] = mypc.m_B_external_particle[2];
+    }
+    else if (mypc.m_B_ext_particle_s=="parse_e_ext_particle_function")
+    {
+        m_type = Parser;
+        m_time = warpx.gett_new(a_pti.GetLevel());
+        m_get_position = GetParticlePosition(a_pti, a_offset);
+        m_xfield_partparser = mypc.m_Bx_particle_parser.get();
+        m_yfield_partparser = mypc.m_By_particle_parser.get();
+        m_zfield_partparser = mypc.m_Bz_particle_parser.get();
+    }
+}
diff --git a/Source/Particles/Gather/Make.package b/Source/Particles/Gather/Make.package
index 19eb98e14..c9a878217 100644
--- a/Source/Particles/Gather/Make.package
+++ b/Source/Particles/Gather/Make.package
@@ -1 +1,3 @@
+CEXE_sources += GetExternalFields.cpp
+
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Particles/Gather
diff --git a/Source/Particles/Gather/ScaleFields.H b/Source/Particles/Gather/ScaleFields.H
new file mode 100644
index 000000000..cc94fac29
--- /dev/null
+++ b/Source/Particles/Gather/ScaleFields.H
@@ -0,0 +1,63 @@
+#ifndef WARPX_PARTICLES_GATHER_SCALEFIELDS_H_
+#define WARPX_PARTICLES_GATHER_SCALEFIELDS_H_
+
+#include "Particles/WarpXParticleContainer.H"
+
+#include <AMReX_REAL.H>
+
+#include <limits>
+
+/** \brief Functor that scales E and B by a factor before pushing the particles.
+ *   This is used for rigid injection.
+ */
+struct ScaleFields
+{
+    amrex::Real m_do_scale;
+    amrex::Real m_dt;
+    amrex::Real m_z_plane_previous;
+    amrex::Real m_vz_ave_boosted;
+    amrex::Real m_v_boost;
+
+    ScaleFields(bool do_scale) noexcept
+        : m_do_scale(do_scale)
+    {}
+
+    ScaleFields (bool do_scale, amrex::Real dt, amrex::Real z_plane_previous,
+                 amrex::Real vz_ave_boosted, amrex::Real v_boost) noexcept
+    : m_do_scale(do_scale), m_dt(dt), m_z_plane_previous(z_plane_previous),
+        m_vz_ave_boosted(vz_ave_boosted), m_v_boost(v_boost)
+    {}
+
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    void operator () (amrex::ParticleReal  /*xp*/,
+                      amrex::ParticleReal  /*yp*/,
+                      amrex::ParticleReal  zp,
+                      amrex::ParticleReal& Exp,
+                      amrex::ParticleReal& Eyp,
+                      amrex::ParticleReal& Ezp,
+                      amrex::ParticleReal& Bxp,
+                      amrex::ParticleReal& Byp,
+                      amrex::ParticleReal& Bzp) const noexcept
+    {
+        using namespace amrex::literals;
+
+        if (not m_do_scale) return;
+
+        // Scale the fields of particles about to cross the injection plane.
+        // This only approximates what should be happening. The particles
+        // should by advanced a fraction of a time step instead.
+        // Scaling the fields is much easier and may be good enough.
+        const amrex::Real dtscale = m_dt - (m_z_plane_previous - zp)/(m_vz_ave_boosted + m_v_boost);
+        if (0._rt < dtscale && dtscale < m_dt)
+        {
+            Exp *= dtscale;
+            Eyp *= dtscale;
+            Ezp *= dtscale;
+            Bxp *= dtscale;
+            Byp *= dtscale;
+            Bzp *= dtscale;
+        }
+    }
+};
+
+#endif
diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H
index 194c28e9e..ecd6717c0 100644
--- a/Source/Particles/MultiParticleContainer.H
+++ b/Source/Particles/MultiParticleContainer.H
@@ -79,15 +79,6 @@ public:
     void InitData ();
 
     ///
-    /// Performs the field gather operation using the input fields E and B, for all the species
-    /// in the MultiParticleContainer. This is the electromagnetic version of the field gather.
-    ///
-    void FieldGather (int lev,
-                      const amrex::MultiFab& Ex, const amrex::MultiFab& Ey,
-                      const amrex::MultiFab& Ez, const amrex::MultiFab& Bx,
-                      const amrex::MultiFab& By, const amrex::MultiFab& Bz);
-
-    ///
     /// This evolves all the particles by one PIC time step, including current deposition, the
     /// field solve, and pushing the particles, for all the species in the MultiParticleContainer.
     /// This is the electromagnetic version.
@@ -125,7 +116,9 @@ public:
     ///
     std::unique_ptr<amrex::MultiFab> GetChargeDensity(int lev, bool local = false);
 
-    void doFieldIonization ();
+    void doFieldIonization (int lev,
+                            const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez,
+                            const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz);
 
     void doCoulombCollisions ();
 
@@ -219,7 +212,13 @@ public:
     /**
     * \brief Performs QED events (Breit-Wheeler process and photon emission)
     */
-    void doQedEvents();
+    void doQedEvents (int lev,
+                      const amrex::MultiFab& Ex,
+                      const amrex::MultiFab& Ey,
+                      const amrex::MultiFab& Ez,
+                      const amrex::MultiFab& Bx,
+                      const amrex::MultiFab& By,
+                      const amrex::MultiFab& Bz);
 #endif
 
     int getSpeciesID (std::string product_str) const;
@@ -230,12 +229,24 @@ protected:
     /**
     * \brief Performs Breit-Wheeler process for the species for which it is enabled
     */
-    void doQedBreitWheeler();
+    void doQedBreitWheeler (int lev,
+                            const amrex::MultiFab& Ex,
+                            const amrex::MultiFab& Ey,
+                            const amrex::MultiFab& Ez,
+                            const amrex::MultiFab& Bx,
+                            const amrex::MultiFab& By,
+                            const amrex::MultiFab& Bz);
 
     /**
     * \brief Performs QED photon emission for the species for which it is enabled
     */
-    void doQedQuantumSync();
+    void doQedQuantumSync (int lev,
+                           const amrex::MultiFab& Ex,
+                           const amrex::MultiFab& Ey,
+                           const amrex::MultiFab& Ez,
+                           const amrex::MultiFab& Bx,
+                           const amrex::MultiFab& By,
+                           const amrex::MultiFab& Bz);
 #endif
 
     // Particle container types
diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index b22a6f79d..fd0980d45 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -313,17 +313,6 @@ MultiParticleContainer::InitData ()
 }
 
 void
-MultiParticleContainer::FieldGather (int lev,
-                                     const MultiFab& Ex, const MultiFab& Ey,
-                                     const MultiFab& Ez, const MultiFab& Bx,
-                                     const MultiFab& By, const MultiFab& Bz)
-{
-    for (auto& pc : allcontainers) {
-        pc->FieldGather(lev, Ex, Ey, Ez, Bx, By, Bz);
-    }
-}
-
-void
 MultiParticleContainer::Evolve (int lev,
                                 const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
                                 const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz,
@@ -634,7 +623,13 @@ MultiParticleContainer::getSpeciesID (std::string product_str) const
 }
 
 void
-MultiParticleContainer::doFieldIonization ()
+MultiParticleContainer::doFieldIonization (int lev,
+                                           const MultiFab& Ex,
+                                           const MultiFab& Ey,
+                                           const MultiFab& Ez,
+                                           const MultiFab& Bx,
+                                           const MultiFab& By,
+                                           const MultiFab& Bz)
 {
     WARPX_PROFILE("MPC::doFieldIonization");
 
@@ -649,31 +644,31 @@ MultiParticleContainer::doFieldIonization ()
         SmartCopyFactory copy_factory(*pc_source, *pc_product);
         auto phys_pc_ptr = static_cast<PhysicalParticleContainer*>(pc_source.get());
 
-        auto Filter    = phys_pc_ptr->getIonizationFunc();
         auto Copy      = copy_factory.getSmartCopy();
         auto Transform = IonizationTransformFunc();
 
         pc_source ->defineAllParticleTiles();
         pc_product->defineAllParticleTiles();
 
-        for (int lev = 0; lev <= pc_source->finestLevel(); ++lev)
-        {
-            const auto info = getMFItInfo(*pc_source, *pc_product);
+        auto info = getMFItInfo(*pc_source, *pc_product);
 
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
-            for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi)
-            {
-                auto& src_tile = pc_source ->ParticlesAt(lev, mfi);
-                auto& dst_tile = pc_product->ParticlesAt(lev, mfi);
+        for (WarpXParIter pti(*pc_source, lev, info); pti.isValid(); ++pti)
+        {
+            auto& src_tile = pc_source ->ParticlesAt(lev, pti);
+            auto& dst_tile = pc_product->ParticlesAt(lev, pti);
 
-                const auto np_dst = dst_tile.numParticles();
-                const auto num_added = filterCopyTransformParticles<1>(dst_tile, src_tile, np_dst,
-                                                                 Filter, Copy, Transform);
+            auto Filter = phys_pc_ptr->getIonizationFunc(pti, lev, Ex.nGrow(),
+                                                         Ex[pti], Ey[pti], Ez[pti],
+                                                         Bx[pti], By[pti], Bz[pti]);
 
-                setNewParticleIDs(dst_tile, np_dst, num_added);
-            }
+            const auto np_dst = dst_tile.numParticles();
+            const auto num_added = filterCopyTransformParticles<1>(dst_tile, src_tile, np_dst,
+                                                                   Filter, Copy, Transform);
+
+            setNewParticleIDs(dst_tile, np_dst, num_added);
         }
     }
 }
@@ -1126,15 +1121,27 @@ MultiParticleContainer::doQEDSchwinger ()
     }
 }
 
-void MultiParticleContainer::doQedEvents()
+void MultiParticleContainer::doQedEvents (int lev,
+                                          const MultiFab& Ex,
+                                          const MultiFab& Ey,
+                                          const MultiFab& Ez,
+                                          const MultiFab& Bx,
+                                          const MultiFab& By,
+                                          const MultiFab& Bz)
 {
     WARPX_PROFILE("MPC::doQedEvents");
 
-    doQedBreitWheeler();
-    doQedQuantumSync();
+    doQedBreitWheeler(lev, Ex, Ey, Ez, Bx, By, Bz);
+    doQedQuantumSync(lev, Ex, Ey, Ez, Bx, By, Bz);
 }
 
-void MultiParticleContainer::doQedBreitWheeler()
+void MultiParticleContainer::doQedBreitWheeler (int lev,
+                                                const MultiFab& Ex,
+                                                const MultiFab& Ey,
+                                                const MultiFab& Ez,
+                                                const MultiFab& Bx,
+                                                const MultiFab& By,
+                                                const MultiFab& Bz)
 {
     WARPX_PROFILE("MPC::doQedBreitWheeler");
 
@@ -1160,40 +1167,48 @@ void MultiParticleContainer::doQedBreitWheeler()
         const auto CopyPos = copy_factory_pos.getSmartCopy();
 
         const auto pair_gen_functor = m_shr_p_bw_engine->build_pair_functor();
-        auto Transform = PairGenerationTransformFunc(pair_gen_functor);
 
         pc_source ->defineAllParticleTiles();
         pc_product_pos->defineAllParticleTiles();
         pc_product_ele->defineAllParticleTiles();
 
-        for (int lev = 0; lev <= pc_source->finestLevel(); ++lev)
-        {
-            const auto info = getMFItInfo(*pc_source, *pc_product_ele, *pc_product_pos);
+        auto info = getMFItInfo(*pc_source, *pc_product_ele, *pc_product_pos);
 
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
-            for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi)
-            {
-                auto& src_tile = pc_source->ParticlesAt(lev, mfi);
-                auto& dst_ele_tile = pc_product_ele->ParticlesAt(lev, mfi);
-                auto& dst_pos_tile = pc_product_pos->ParticlesAt(lev, mfi);
-
-                const auto np_dst_ele = dst_ele_tile.numParticles();
-                const auto np_dst_pos = dst_pos_tile.numParticles();
-                const auto num_added = filterCopyTransformParticles<1>(
-                    dst_ele_tile, dst_pos_tile,
-                    src_tile, np_dst_ele, np_dst_pos,
-                    Filter, CopyEle, CopyPos, Transform);
-
-                setNewParticleIDs(dst_ele_tile, np_dst_ele, num_added);
-                setNewParticleIDs(dst_pos_tile, np_dst_pos, num_added);
-            }
+        for (WarpXParIter pti(*pc_source, lev, info); pti.isValid(); ++pti)
+        {
+            auto Transform = PairGenerationTransformFunc(pair_gen_functor,
+                                                         pti, lev, Ex.nGrow(),
+                                                         Ex[pti], Ey[pti], Ez[pti],
+                                                         Bx[pti], By[pti], Bz[pti],
+                                                         pc_source->get_v_galilean());
+
+            auto& src_tile = pc_source->ParticlesAt(lev, pti);
+            auto& dst_ele_tile = pc_product_ele->ParticlesAt(lev, pti);
+            auto& dst_pos_tile = pc_product_pos->ParticlesAt(lev, pti);
+
+            const auto np_dst_ele = dst_ele_tile.numParticles();
+            const auto np_dst_pos = dst_pos_tile.numParticles();
+            const auto num_added = filterCopyTransformParticles<1>(
+                                                      dst_ele_tile, dst_pos_tile,
+                                                      src_tile, np_dst_ele, np_dst_pos,
+                                                      Filter, CopyEle, CopyPos, Transform);
+
+            setNewParticleIDs(dst_ele_tile, np_dst_ele, num_added);
+            setNewParticleIDs(dst_pos_tile, np_dst_pos, num_added);
         }
-   }
+    }
 }
 
-void MultiParticleContainer::doQedQuantumSync()
+void MultiParticleContainer::doQedQuantumSync (int lev,
+                                               const MultiFab& Ex,
+                                               const MultiFab& Ey,
+                                               const MultiFab& Ez,
+                                               const MultiFab& Bx,
+                                               const MultiFab& By,
+                                               const MultiFab& Bz)
 {
     WARPX_PROFILE("MPC::doQedEvents::doQedQuantumSync");
 
@@ -1204,7 +1219,7 @@ void MultiParticleContainer::doQedQuantumSync()
     for (auto& pc_source : allcontainers){
         if(!pc_source->has_quantum_sync()){ continue; }
 
-                // Get product species
+        // Get product species
         auto& pc_product_phot =
             allcontainers[pc_source->m_qed_quantum_sync_phot_product];
 
@@ -1215,42 +1230,41 @@ void MultiParticleContainer::doQedQuantumSync()
         const auto Filter   = phys_pc_ptr->getPhotonEmissionFilterFunc();
         const auto CopyPhot = copy_factory_phot.getSmartCopy();
 
-        auto Transform = PhotonEmissionTransformFunc(
-            m_shr_p_qs_engine->build_optical_depth_functor(),
-            pc_source->particle_runtime_comps["optical_depth_QSR"],
-             m_shr_p_qs_engine->build_phot_em_functor());
-
         pc_source ->defineAllParticleTiles();
         pc_product_phot->defineAllParticleTiles();
 
-        for (int lev = 0; lev <= pc_source->finestLevel(); ++lev)
-        {
-            const auto info = getMFItInfo(*pc_source, *pc_product_phot);
+        auto info = getMFItInfo(*pc_source, *pc_product_phot);
 
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
-            for (MFIter mfi = pc_source->MakeMFIter(lev, info); mfi.isValid(); ++mfi)
-            {
-                auto& src_tile = pc_source->ParticlesAt(lev, mfi);
-                auto& dst_tile = pc_product_phot->ParticlesAt(lev, mfi);
+        for (WarpXParIter pti(*pc_source, lev, info); pti.isValid(); ++pti)
+        {
+            auto Transform = PhotonEmissionTransformFunc(
+                  m_shr_p_qs_engine->build_optical_depth_functor(),
+                  pc_source->particle_runtime_comps["optical_depth_QSR"],
+                  m_shr_p_qs_engine->build_phot_em_functor(),
+                  pti, lev, Ex.nGrow(),
+                  Ex[pti], Ey[pti], Ez[pti],
+                  Bx[pti], By[pti], Bz[pti],
+                  pc_source->get_v_galilean());
 
-                const auto np_dst = dst_tile.numParticles();
+            auto& src_tile = pc_source->ParticlesAt(lev, pti);
+            auto& dst_tile = pc_product_phot->ParticlesAt(lev, pti);
 
-                const auto num_added =
-                    filterCopyTransformParticles<1>(dst_tile, src_tile, np_dst,
-                        Filter, CopyPhot, Transform);
+            const auto np_dst = dst_tile.numParticles();
 
-                setNewParticleIDs(dst_tile, np_dst, num_added);
+            const auto num_added =
+                filterCopyTransformParticles<1>(dst_tile, src_tile, np_dst,
+                                                Filter, CopyPhot, Transform);
 
-                cleanLowEnergyPhotons(
-                    dst_tile, np_dst, num_added,
-                    m_quantum_sync_photon_creation_energy_threshold);
+            setNewParticleIDs(dst_tile, np_dst, num_added);
 
-            }
+            cleanLowEnergyPhotons(
+                                  dst_tile, np_dst, num_added,
+                                  m_quantum_sync_photon_creation_energy_threshold);
         }
     }
-
 }
 
 void MultiParticleContainer::CheckQEDProductSpecies()
diff --git a/Source/Particles/ParticleCreation/DefaultInitialization.H b/Source/Particles/ParticleCreation/DefaultInitialization.H
index 9a8d9f7ee..0278c7118 100644
--- a/Source/Particles/ParticleCreation/DefaultInitialization.H
+++ b/Source/Particles/ParticleCreation/DefaultInitialization.H
@@ -41,12 +41,6 @@ static std::map<std::string, InitializationPolicy> initialization_policies = {
     {"ux",    InitializationPolicy::Zero },
     {"uy",    InitializationPolicy::Zero },
     {"uz",    InitializationPolicy::Zero },
-    {"Ex",    InitializationPolicy::Zero },
-    {"Ey",    InitializationPolicy::Zero },
-    {"Ez",    InitializationPolicy::Zero },
-    {"Bx",    InitializationPolicy::Zero },
-    {"By",    InitializationPolicy::Zero },
-    {"Bz",    InitializationPolicy::Zero },
 #ifdef WARPX_DIM_RZ
     {"theta", InitializationPolicy::Zero},
 #endif
diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H
index 4cdffd385..ad89ca0a5 100644
--- a/Source/Particles/PhotonParticleContainer.H
+++ b/Source/Particles/PhotonParticleContainer.H
@@ -55,7 +55,18 @@ public:
                          DtType a_dt_type=DtType::Full) override;
 
     virtual void PushPX(WarpXParIter& pti,
-                        amrex::Real dt, DtType a_dt_type=DtType::Full) override;
+                        amrex::FArrayBox const * exfab,
+                        amrex::FArrayBox const * eyfab,
+                        amrex::FArrayBox const * ezfab,
+                        amrex::FArrayBox const * bxfab,
+                        amrex::FArrayBox const * byfab,
+                        amrex::FArrayBox const * bzfab,
+                        const int ngE, const int /*e_is_nodal*/,
+                        const long offset,
+                        const long np_to_push,
+                        int lev, int gather_lev,
+                        amrex::Real dt, ScaleFields scaleFields,
+                        DtType a_dt_type) override;
 
     // Do nothing
     virtual void PushP (int lev,
@@ -84,19 +95,6 @@ public:
                                 int lev,
                                 int depos_lev,
                                 amrex::Real dt) override {};
-
-#ifdef WARPX_QED
-    /**
-     * This function evolves the optical depth of the photons if QED effects
-     * are enabled.
-     * @param[in,out] pti particle iterator (optical depth will be modified)
-     * @param[in] dt temporal step
-     */
-    virtual void EvolveOpticalDepth(WarpXParIter& pti,
-        amrex::Real dt) override;
-
-#endif
-
 };
 
 #endif // #ifndef WARPX_PhotonParticleContainer_H_
diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp
index bd80b8e43..abf56dd7c 100644
--- a/Source/Particles/PhotonParticleContainer.cpp
+++ b/Source/Particles/PhotonParticleContainer.cpp
@@ -13,9 +13,11 @@
 #include "Particles/Pusher/UpdatePositionPhoton.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/Pusher/CopyParticleAttribs.H"
+#include "Particles/Gather/FieldGather.H"
+#include "Particles/Gather/GetExternalFields.H"
 
 #ifdef _OPENMP
-#   include <omp.h>
+#include <omp.h>
 #endif
 
 #include <limits>
@@ -63,10 +65,35 @@ void PhotonParticleContainer::InitData()
 }
 
 void
-PhotonParticleContainer::PushPX(WarpXParIter& pti, Real dt, DtType a_dt_type)
+PhotonParticleContainer::PushPX (WarpXParIter& pti,
+                                 amrex::FArrayBox const * exfab,
+                                 amrex::FArrayBox const * eyfab,
+                                 amrex::FArrayBox const * ezfab,
+                                 amrex::FArrayBox const * bxfab,
+                                 amrex::FArrayBox const * byfab,
+                                 amrex::FArrayBox const * bzfab,
+                                 const int ngE, const int /*e_is_nodal*/,
+                                 const long offset,
+                                 const long np_to_push,
+                                 int lev, int gather_lev,
+                                 amrex::Real dt, ScaleFields /*scaleFields*/, DtType a_dt_type)
 {
+    // Get cell size on gather_lev
+    const std::array<Real,3>& dx = WarpX::CellSize(std::max(gather_lev,0));
+
+    // Get box from which field is gathered.
+    // If not gathering from the finest level, the box is coarsened.
+    amrex::Box box;
+    if (lev == gather_lev) {
+        box = pti.tilebox();
+    } else {
+        const IntVect& ref_ratio = WarpX::RefRatio(gather_lev);
+        box = amrex::coarsen(pti.tilebox(),ref_ratio);
+    }
+
+    // Add guard cells to the box.
+    box.grow(ngE);
 
-    // This wraps the momentum and position advance so that inheritors can modify the call.
     auto& attribs = pti.GetAttribs();
 
     // Extract pointers to the different particle quantities
@@ -74,19 +101,91 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, Real dt, DtType a_dt_type)
     ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
     ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
 
+#ifdef WARPX_QED
+    AMREX_ASSERT(has_breit_wheeler());
+
+    BreitWheelerEvolveOpticalDepth evolve_opt;
+    amrex::Real* AMREX_RESTRICT p_optical_depth_BW = nullptr;
+    const bool local_has_breit_wheeler = has_breit_wheeler();
+    if (local_has_breit_wheeler) {
+        evolve_opt = m_shr_p_bw_engine->build_evolve_functor();
+        p_optical_depth_BW = pti.GetAttribs(particle_comps["optical_depth_BW"]).dataPtr();
+    }
+
+    const auto me = PhysConst::m_e;
+#endif
+
     auto copyAttribs = CopyParticleAttribs(pti, tmp_particle_data);
     int do_copy = (WarpX::do_back_transformed_diagnostics &&
                    do_back_transformed_diagnostics && a_dt_type!=DtType::SecondHalf);
 
-    const auto GetPosition = GetParticlePosition(pti);
-          auto SetPosition = SetParticlePosition(pti);
+    const auto GetPosition = GetParticlePosition(pti, offset);
+    auto SetPosition = SetParticlePosition(pti, offset);
+
+    const auto getExternalE = GetExternalEField(pti, offset);
+    const auto getExternalB = GetExternalBField(pti, offset);
+
+    // Lower corner of tile box physical domain (take into account Galilean shift)
+    amrex::Real cur_time = WarpX::GetInstance().gett_new(lev);
+    const auto& time_of_last_gal_shift = WarpX::GetInstance().time_of_last_gal_shift;
+    amrex::Real time_shift = (cur_time - time_of_last_gal_shift);
+    amrex::Array<amrex::Real,3> galilean_shift = { v_galilean[0]*time_shift, v_galilean[1]*time_shift, v_galilean[2]*time_shift };
+    const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, galilean_shift, gather_lev);
+
+    const Dim3 lo = lbound(box);
+
+    int l_lower_order_in_v = WarpX::l_lower_order_in_v;
+    int nox = WarpX::nox;
+    int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+    amrex::GpuArray<amrex::Real, 3> dx_arr = {dx[0], dx[1], dx[2]};
+    amrex::GpuArray<amrex::Real, 3> xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+    amrex::Array4<const amrex::Real> const& ex_arr = exfab->array();
+    amrex::Array4<const amrex::Real> const& ey_arr = eyfab->array();
+    amrex::Array4<const amrex::Real> const& ez_arr = ezfab->array();
+    amrex::Array4<const amrex::Real> const& bx_arr = bxfab->array();
+    amrex::Array4<const amrex::Real> const& by_arr = byfab->array();
+    amrex::Array4<const amrex::Real> const& bz_arr = bzfab->array();
+
+    amrex::IndexType const ex_type = exfab->box().ixType();
+    amrex::IndexType const ey_type = eyfab->box().ixType();
+    amrex::IndexType const ez_type = ezfab->box().ixType();
+    amrex::IndexType const bx_type = bxfab->box().ixType();
+    amrex::IndexType const by_type = byfab->box().ixType();
+    amrex::IndexType const bz_type = bzfab->box().ixType();
 
     amrex::ParallelFor(
-        pti.numParticles(),
+        np_to_push,
         [=] AMREX_GPU_DEVICE (long i) {
             if (do_copy) copyAttribs(i);
             ParticleReal x, y, z;
             GetPosition(i, x, y, z);
+
+            amrex::ParticleReal Exp, Eyp, Ezp;
+            getExternalE(i, Exp, Eyp, Ezp);
+
+            amrex::ParticleReal Bxp, Byp, Bzp;
+            getExternalB(i, Bxp, Byp, Bzp);
+
+            // first gather E and B to the particle positions
+            doGatherShapeN(x, y, z, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                           ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                           ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                           dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes,
+                           nox, l_lower_order_in_v);
+
+#ifdef WARPX_QED
+            if (local_has_breit_wheeler) {
+                const ParticleReal px = me * ux[i];
+                const ParticleReal py = me * uy[i];
+                const ParticleReal pz = me * uz[i];
+
+                bool has_event_happened = evolve_opt(px, py, pz, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                                                     dt, p_optical_depth_BW[i]);
+            }
+#endif
+
             UpdatePositionPhoton( x, y, z, ux[i], uy[i], uz[i], dt );
             SetPosition(i, x, y, z);
         }
@@ -118,49 +217,3 @@ PhotonParticleContainer::Evolve (int lev,
                                        t, dt);
 
 }
-
-#ifdef WARPX_QED
-
-void
-PhotonParticleContainer::EvolveOpticalDepth(
-    WarpXParIter& pti,amrex::Real dt)
-{
-    if(!has_breit_wheeler())
-        return;
-
-    auto& attribs = pti.GetAttribs();
-    ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
-    ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
-    ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr();
-
-    BreitWheelerEvolveOpticalDepth evolve_opt =
-        m_shr_p_bw_engine->build_evolve_functor();
-
-    amrex::Real* AMREX_RESTRICT p_optical_depth_BW =
-        pti.GetAttribs(particle_comps["optical_depth_BW"]).dataPtr();
-
-    const auto me = PhysConst::m_e;
-
-    amrex::ParallelFor(
-        pti.numParticles(),
-        [=] AMREX_GPU_DEVICE (long i) {
-            const ParticleReal px = me * ux[i];
-            const ParticleReal py = me * uy[i];
-            const ParticleReal pz = me * uz[i];
-
-            bool has_event_happened = evolve_opt(
-                px, py, pz,
-                Ex[i], Ey[i], Ez[i],
-                Bx[i], By[i], Bz[i],
-                dt, p_optical_depth_BW[i]);
-        }
-        );
-}
-
-#endif
diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H
index afa7a3b85..5d8102d09 100644
--- a/Source/Particles/PhysicalParticleContainer.H
+++ b/Source/Particles/PhysicalParticleContainer.H
@@ -14,6 +14,7 @@
 #include "WarpXParticleContainer.H"
 #include "Filter/NCIGodfreyFilter.H"
 #include "Particles/ElementaryProcess/Ionization.H"
+#include "Particles/Gather/ScaleFields.H"
 
 #ifdef WARPX_QED
 #    include "Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H"
@@ -63,33 +64,6 @@ public:
 
     void InitIonizationModule ();
 
-    virtual void FieldGather (int lev,
-                              const amrex::MultiFab& Ex,
-                              const amrex::MultiFab& Ey,
-                              const amrex::MultiFab& Ez,
-                              const amrex::MultiFab& Bx,
-                              const amrex::MultiFab& By,
-                              const amrex::MultiFab& Bz) final;
-
-    void FieldGather (WarpXParIter& pti,
-                      RealVector& Exp,
-                      RealVector& Eyp,
-                      RealVector& Ezp,
-                      RealVector& Bxp,
-                      RealVector& Byp,
-                      RealVector& Bzp,
-                      amrex::FArrayBox const * exfab,
-                      amrex::FArrayBox const * eyfab,
-                      amrex::FArrayBox const * ezfab,
-                      amrex::FArrayBox const * bxfab,
-                      amrex::FArrayBox const * byfab,
-                      amrex::FArrayBox const * bzfab,
-                      const int ngE, const int e_is_nodal,
-                      const long offset,
-                      const long np_to_gather,
-                      int lev,
-                      int depos_lev);
-
     /**
      * \brief Evolve is the central function PhysicalParticleContainer that
      * advances plasma particles for a time dt (typically one timestep).
@@ -148,7 +122,19 @@ public:
                          amrex::Real dt,
                          DtType a_dt_type=DtType::Full) override;
 
-    virtual void PushPX (WarpXParIter& pti, amrex::Real dt, DtType a_dt_type=DtType::Full);
+    virtual void PushPX (WarpXParIter& pti,
+                         amrex::FArrayBox const * exfab,
+                         amrex::FArrayBox const * eyfab,
+                         amrex::FArrayBox const * ezfab,
+                         amrex::FArrayBox const * bxfab,
+                         amrex::FArrayBox const * byfab,
+                         amrex::FArrayBox const * bzfab,
+                         const int ngE, const int /*e_is_nodal*/,
+                         const long offset,
+                         const long np_to_push,
+                         int lev, int gather_lev,
+                         amrex::Real dt, ScaleFields scaleFields,
+                         DtType a_dt_type=DtType::Full);
 
     virtual void PushP (int lev, amrex::Real dt,
                         const amrex::MultiFab& Ex,
@@ -175,7 +161,15 @@ public:
 
     void SplitParticles (int lev);
 
-    IonizationFilterFunc getIonizationFunc ();
+    IonizationFilterFunc getIonizationFunc (const WarpXParIter& pti,
+                                            int lev,
+                                            int ngE,
+                                            const amrex::FArrayBox& Ex,
+                                            const amrex::FArrayBox& Ey,
+                                            const amrex::FArrayBox& Ez,
+                                            const amrex::FArrayBox& Bx,
+                                            const amrex::FArrayBox& By,
+                                            const amrex::FArrayBox& Bz);
 
     // Inject particles in Box 'part_box'
     virtual void AddParticles (int lev);
@@ -304,15 +298,6 @@ public:
         (std::shared_ptr<QuantumSynchrotronEngine> ptr) override;
     //__________
 
-    /**
-     * This function evolves the optical depth of the particles if QED effects
-     * are enabled.
-     * @param[in,out] pti particle iterator (optical depth will be modified)
-     * @param[in] dt temporal step
-     */
-    virtual void EvolveOpticalDepth (WarpXParIter& pti,
-        amrex::Real dt);
-
     PhotonEmissionFilterFunc getPhotonEmissionFilterFunc ();
 
     PairGenerationFilterFunc getPairGenerationFilterFunc ();
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index a0016d436..07c307c83 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -899,69 +899,6 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
 }
 
 void
-PhysicalParticleContainer::FieldGather (int lev,
-                                        const amrex::MultiFab& Ex,
-                                        const amrex::MultiFab& Ey,
-                                        const amrex::MultiFab& Ez,
-                                        const amrex::MultiFab& Bx,
-                                        const amrex::MultiFab& By,
-                                        const amrex::MultiFab& Bz)
-{
-    BL_ASSERT(OnSameGrids(lev,Ex));
-
-    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
-
-#ifdef _OPENMP
-#pragma omp parallel
-#endif
-    {
-        for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
-        {
-            if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
-            {
-                amrex::Gpu::synchronize();
-            }
-            Real wt = amrex::second();
-
-            auto& attribs = pti.GetAttribs();
-
-            auto& Exp = attribs[PIdx::Ex];
-            auto& Eyp = attribs[PIdx::Ey];
-            auto& Ezp = attribs[PIdx::Ez];
-            auto& Bxp = attribs[PIdx::Bx];
-            auto& Byp = attribs[PIdx::By];
-            auto& Bzp = attribs[PIdx::Bz];
-
-            const long np = pti.numParticles();
-
-            // Data on the grid
-            const FArrayBox& exfab = Ex[pti];
-            const FArrayBox& eyfab = Ey[pti];
-            const FArrayBox& ezfab = Ez[pti];
-            const FArrayBox& bxfab = Bx[pti];
-            const FArrayBox& byfab = By[pti];
-            const FArrayBox& bzfab = Bz[pti];
-
-            //
-            // Field Gather
-            //
-            int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
-            FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
-                        Ex.nGrow(), e_is_nodal,
-                        0, np, lev, lev);
-
-            if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
-            {
-                amrex::Gpu::synchronize();
-                wt = amrex::second() - wt;
-                amrex::HostDevice::Atomic::Add( &(*cost)[pti.index()], wt);
-            }
-        }
-    }
-}
-
-void
 PhysicalParticleContainer::Evolve (int lev,
                                    const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
                                    const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz,
@@ -973,10 +910,7 @@ PhysicalParticleContainer::Evolve (int lev,
                                    Real /*t*/, Real dt, DtType a_dt_type)
 {
     WARPX_PROFILE("PPC::Evolve()");
-    WARPX_PROFILE_VAR_NS("PPC::Evolve::Copy", blp_copy);
-    WARPX_PROFILE_VAR_NS("PPC::FieldGather", blp_fg);
-    WARPX_PROFILE_VAR_NS("PPC::EvolveOpticalDepth", blp_ppc_qed_ev);
-    WARPX_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp);
+    WARPX_PROFILE_VAR_NS("PPC::GatherAndPush", blp_fg);
 
     BL_ASSERT(OnSameGrids(lev,jx));
 
@@ -1029,12 +963,6 @@ PhysicalParticleContainer::Evolve (int lev,
             auto& uxp = attribs[PIdx::ux];
             auto& uyp = attribs[PIdx::uy];
             auto& uzp = attribs[PIdx::uz];
-            auto& Exp = attribs[PIdx::Ex];
-            auto& Eyp = attribs[PIdx::Ey];
-            auto& Ezp = attribs[PIdx::Ez];
-            auto& Bxp = attribs[PIdx::Bx];
-            auto& Byp = attribs[PIdx::By];
-            auto& Bzp = attribs[PIdx::Bz];
 
             const long np = pti.numParticles();
 
@@ -1103,13 +1031,13 @@ PhysicalParticleContainer::Evolve (int lev,
                 int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
 
                 //
-                // Field Gather of Aux Data (i.e., the full solution)
+                // Gather and push for particles not in the buffer
                 //
                 WARPX_PROFILE_VAR_START(blp_fg);
-                FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                            exfab, eyfab, ezfab, bxfab, byfab, bzfab,
-                            Ex.nGrow(), e_is_nodal,
-                            0, np_gather, lev, lev);
+                PushPX(pti, exfab, eyfab, ezfab,
+                       bxfab, byfab, bzfab,
+                       Ex.nGrow(), e_is_nodal,
+                       0, np_gather, lev, lev, dt, ScaleFields(false), a_dt_type);
 
                 if (np_gather < np)
                 {
@@ -1138,34 +1066,17 @@ PhysicalParticleContainer::Evolve (int lev,
                                        cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab);
                     }
 
-                    // Field gather for particles in gather buffers
+                    // Field gather and push for particles in gather buffers
                     e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal();
-                    FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                                cexfab, ceyfab, cezfab,
-                                cbxfab, cbyfab, cbzfab,
-                                cEx->nGrow(), e_is_nodal,
-                                nfine_gather, np-nfine_gather,
-                                lev, lev-1);
+                    PushPX(pti, cexfab, ceyfab, cezfab,
+                           cbxfab, cbyfab, cbzfab,
+                           cEx->nGrow(), e_is_nodal,
+                           nfine_gather, np-nfine_gather,
+                           lev, lev-1, dt, ScaleFields(false), a_dt_type);
                 }
 
                 WARPX_PROFILE_VAR_STOP(blp_fg);
 
-#ifdef WARPX_QED
-                //
-                //Evolve Optical Depth
-                //
-                WARPX_PROFILE_VAR_START(blp_ppc_qed_ev);
-                EvolveOpticalDepth(pti, dt);
-                WARPX_PROFILE_VAR_STOP(blp_ppc_qed_ev);
-#endif
-
-                //
-                // Particle Push
-                //
-                WARPX_PROFILE_VAR_START(blp_ppc_pp);
-                PushPX(pti, dt, a_dt_type);
-                WARPX_PROFILE_VAR_STOP(blp_ppc_pp);
-
                 //
                 // Current Deposition (only needed for electromagnetic solver)
                 //
@@ -1466,131 +1377,24 @@ PhysicalParticleContainer::SplitParticles (int lev)
 }
 
 void
-PhysicalParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type)
-{
-
-    // This wraps the momentum and position advance so that inheritors can modify the call.
-    auto& attribs = pti.GetAttribs();
-    // Extract pointers to the different particle quantities
-
-    const auto GetPosition = GetParticlePosition(pti);
-          auto SetPosition = SetParticlePosition(pti);
-
-    ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
-    ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
-    ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr();
-
-    auto copyAttribs = CopyParticleAttribs(pti, tmp_particle_data);
-    int do_copy = (WarpX::do_back_transformed_diagnostics &&
-                          do_back_transformed_diagnostics &&
-                   (a_dt_type!=DtType::SecondHalf));
-
-    int* AMREX_RESTRICT ion_lev = nullptr;
-    if (do_field_ionization){
-        ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
-    }
-
-    // Loop over the particles and update their momentum
-    const Real q = this->charge;
-    const Real m = this-> mass;
-
-    const auto pusher_algo = WarpX::particle_pusher_algo;
-    const auto do_crr = do_classical_radiation_reaction;
-#ifdef WARPX_QED
-    const auto do_sync = m_do_qed_quantum_sync;
-    amrex::Real t_chi_max = 0.0;
-    if (do_sync) t_chi_max = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min;
-#endif
-
-    amrex::ParallelFor(pti.numParticles(),
-                       [=] AMREX_GPU_DEVICE (long i) {
-                           doParticlePush(GetPosition, SetPosition, copyAttribs, i,
-                                          ux[i], uy[i], uz[i],
-                                          Ex[i], Ey[i], Ez[i],
-                                          Bx[i], By[i], Bz[i],
-                                          ion_lev ? ion_lev[i] : 0,
-                                          m, q, pusher_algo, do_crr, do_copy,
-#ifdef WARPX_QED
-                                          do_sync,
-                                          t_chi_max,
-#endif
-                                          dt);
-                       });
-}
-
-#ifdef WARPX_QED
-void PhysicalParticleContainer::EvolveOpticalDepth(
-    WarpXParIter& pti, amrex::Real dt)
-{
-    if(!has_quantum_sync())
-        return;
-
-    QuantumSynchrotronEvolveOpticalDepth evolve_opt =
-        m_shr_p_qs_engine->build_evolve_functor();
-
-    auto& attribs = pti.GetAttribs();
-    const ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Ex = attribs[PIdx::Ex].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Ey = attribs[PIdx::Ey].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Ez = attribs[PIdx::Ez].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Bx = attribs[PIdx::Bx].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr();
-    const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr();
-
-    ParticleReal* const AMREX_RESTRICT p_optical_depth_QSR =
-        pti.GetAttribs(particle_comps["optical_depth_QSR"]).dataPtr();
-
-    const ParticleReal m = this->mass;
-
-    amrex::ParallelFor(pti.numParticles(),
-            [=] AMREX_GPU_DEVICE (long i) {
-                const ParticleReal px = m * ux[i];
-                const ParticleReal py = m * uy[i];
-                const ParticleReal pz = m * uz[i];
-
-                bool has_event_happened = evolve_opt(
-                    px, py, pz,
-                    Ex[i], Ey[i], Ez[i],
-                    Bx[i], By[i], Bz[i],
-                    dt, p_optical_depth_QSR[i]);
-            }
-    );
-
-}
-#endif
-
-void
-PhysicalParticleContainer::PushP (
-    int lev, Real dt,
-    const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
-    const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz)
+PhysicalParticleContainer::PushP (int lev, Real dt,
+                                  const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez,
+                                  const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz)
 {
     WARPX_PROFILE("PhysicalParticleContainer::PushP");
 
     if (do_not_push) return;
 
+    const std::array<amrex::Real,3>& dx = WarpX::CellSize(std::max(lev,0));
+
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
     {
         for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
         {
-            auto& attribs = pti.GetAttribs();
-
-            auto& Exp = attribs[PIdx::Ex];
-            auto& Eyp = attribs[PIdx::Ey];
-            auto& Ezp = attribs[PIdx::Ez];
-            auto& Bxp = attribs[PIdx::Bx];
-            auto& Byp = attribs[PIdx::By];
-            auto& Bzp = attribs[PIdx::Bz];
+            amrex::Box box = pti.tilebox();
+            box.grow(Ex.nGrow());
 
             const long np = pti.numParticles();
 
@@ -1602,83 +1406,100 @@ PhysicalParticleContainer::PushP (
             const FArrayBox& byfab = By[pti];
             const FArrayBox& bzfab = Bz[pti];
 
-            int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
-            FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
-                        Ex.nGrow(), e_is_nodal,
-                        0, np, lev, lev);
+            const auto getPosition = GetParticlePosition(pti);
+                  auto setPosition = SetParticlePosition(pti);
+
+            const auto getExternalE = GetExternalEField(pti);
+            const auto getExternalB = GetExternalBField(pti);
+
+            const auto& xyzmin = WarpX::GetInstance().LowerCornerWithGalilean(box,v_galilean,lev);
+
+            const Dim3 lo = lbound(box);
 
-            // This wraps the momentum advance so that inheritors can modify the call.
-            // Extract pointers to the different particle quantities
+            int l_lower_order_in_v = WarpX::l_lower_order_in_v;
+            int nox = WarpX::nox;
+            int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+            amrex::GpuArray<amrex::Real, 3> dx_arr = {dx[0], dx[1], dx[2]};
+            amrex::GpuArray<amrex::Real, 3> xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+            amrex::Array4<const amrex::Real> const& ex_arr = exfab.array();
+            amrex::Array4<const amrex::Real> const& ey_arr = eyfab.array();
+            amrex::Array4<const amrex::Real> const& ez_arr = ezfab.array();
+            amrex::Array4<const amrex::Real> const& bx_arr = bxfab.array();
+            amrex::Array4<const amrex::Real> const& by_arr = byfab.array();
+            amrex::Array4<const amrex::Real> const& bz_arr = bzfab.array();
+
+            amrex::IndexType const ex_type = exfab.box().ixType();
+            amrex::IndexType const ey_type = eyfab.box().ixType();
+            amrex::IndexType const ez_type = ezfab.box().ixType();
+            amrex::IndexType const bx_type = bxfab.box().ixType();
+            amrex::IndexType const by_type = byfab.box().ixType();
+            amrex::IndexType const bz_type = bzfab.box().ixType();
+
+            auto& attribs = pti.GetAttribs();
             ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
             ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
             ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Expp = Exp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Eypp = Eyp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Ezpp = Ezp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Bxpp = Bxp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Bypp = Byp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr();
-
-            // Loop over the particles and update their momentum
-            const Real q = this->charge;
-            const Real m = this-> mass;
 
             int* AMREX_RESTRICT ion_lev = nullptr;
-            if (do_field_ionization){
+            if (do_field_ionization) {
                 ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
             }
 
-            //Assumes that all consistency checks have been done at initialization
-            if(do_classical_radiation_reaction){
-                 amrex::ParallelFor(pti.numParticles(),
-                    [=] AMREX_GPU_DEVICE (long i){
-                        Real qp = q;
-                        if (ion_lev){ qp *= ion_lev[i]; }
-                        UpdateMomentumBorisWithRadiationReaction(
-                            ux[i], uy[i], uz[i],
-                            Expp[i], Eypp[i], Ezpp[i],
-                            Bxpp[i], Bypp[i], Bzpp[i],
-                            qp, m, dt);
-                    }
-                );
-            } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
-                amrex::ParallelFor(pti.numParticles(),
-                    [=] AMREX_GPU_DEVICE (long i) {
-                        Real qp = q;
-                        if (ion_lev){ qp *= ion_lev[i]; }
-                        UpdateMomentumBoris(
-                            ux[i], uy[i], uz[i],
-                            Expp[i], Eypp[i], Ezpp[i],
-                            Bxpp[i], Bypp[i], Bzpp[i],
-                            qp, m, dt);
-                    }
-                );
-            } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay){
-                amrex::ParallelFor(pti.numParticles(),
-                    [=] AMREX_GPU_DEVICE (long i) {
-                        Real qp = q;
-                        if (ion_lev){ qp *= ion_lev[i]; }
-                        UpdateMomentumVay(
-                            ux[i], uy[i], uz[i],
-                            Expp[i], Eypp[i], Ezpp[i],
-                            Bxpp[i], Bypp[i], Bzpp[i],
-                            qp, m, dt);
-                    }
-                );
-            } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary){
-                amrex::ParallelFor(pti.numParticles(),
-                    [=] AMREX_GPU_DEVICE (long i) {
-                        UpdateMomentumHigueraCary( ux[i], uy[i], uz[i],
-                            Expp[i], Eypp[i], Ezpp[i],
-                            Bxpp[i], Bypp[i], Bzpp[i],
-                            q, m, dt);
-                    }
-                );
-            } else {
-              amrex::Abort("Unknown particle pusher");
-            }
+            // Loop over the particles and update their momentum
+            const amrex::Real q = this->charge;
+            const amrex::Real m = this-> mass;
+
+            const auto pusher_algo = WarpX::particle_pusher_algo;
+            const auto do_crr = do_classical_radiation_reaction;
 
+            amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip)
+            {
+                amrex::ParticleReal xp, yp, zp;
+                getPosition(ip, xp, yp, zp);
+
+                amrex::ParticleReal Exp = 0._rt, Eyp = 0._rt, Ezp = 0._rt;
+                getExternalE(ip, Exp, Eyp, Ezp);
+
+                amrex::ParticleReal Bxp = 0._rt, Byp = 0._rt, Bzp = 0._rt;
+                getExternalB(ip, Bxp, Byp, Bzp);
+
+                // first gather E and B to the particle positions
+                doGatherShapeN(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                               ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                               ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                               dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes,
+                               nox, l_lower_order_in_v);
+
+                if (do_crr) {
+                    amrex::Real qp = q;
+                    if (ion_lev) { qp *= ion_lev[ip]; }
+                    UpdateMomentumBorisWithRadiationReaction(ux[ip], uy[ip], uz[ip],
+                                                             Exp, Eyp, Ezp, Bxp,
+                                                             Byp, Bzp, qp, m, dt);
+                } else if (pusher_algo == ParticlePusherAlgo::Boris) {
+                    amrex::Real qp = q;
+                    if (ion_lev) { qp *= ion_lev[ip]; }
+                    UpdateMomentumBoris( ux[ip], uy[ip], uz[ip],
+                                         Exp, Eyp, Ezp, Bxp,
+                                         Byp, Bzp, qp, m, dt);
+                } else if (pusher_algo == ParticlePusherAlgo::Vay) {
+                    amrex::Real qp = q;
+                    if (ion_lev){ qp *= ion_lev[ip]; }
+                    UpdateMomentumVay( ux[ip], uy[ip], uz[ip],
+                                       Exp, Eyp, Ezp, Bxp,
+                                       Byp, Bzp, qp, m, dt);
+                } else if (pusher_algo == ParticlePusherAlgo::HigueraCary) {
+                    amrex::Real qp = q;
+                    if (ion_lev){ qp *= ion_lev[ip]; }
+                    UpdateMomentumHigueraCary( ux[ip], uy[ip], uz[ip],
+                                               Exp, Eyp, Ezp, Bxp,
+                                               Byp, Bzp, qp, m, dt);
+                } else {
+                    amrex::Abort("Unknown particle pusher");
+                }
+            });
         }
     }
 }
@@ -1897,44 +1718,30 @@ PhysicalParticleContainer::ContinuousInjection (const RealBox& injection_box)
     AddPlasma(lev, injection_box);
 }
 
-/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab,
- * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp.
- * \param Exp-Bzp: fields on particles.
- * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti
- * \param ngE: number of guard cells for E
- * \param e_is_nodal: 0 if E is staggered, 1 if E is nodal
- * \param offset: index of first particle for which fields are gathered
- * \param np_to_gather: number of particles onto which fields are gathered
- * \param lev: level on which particles are located
- * \param gather_lev: level from which particles gather fields (lev-1) for
-          particles in buffers.
+/* \brief Perform the field gather and particle push operations in one fused kernel
+ *
  */
 void
-PhysicalParticleContainer::FieldGather (WarpXParIter& pti,
-                                        RealVector& Exp,
-                                        RealVector& Eyp,
-                                        RealVector& Ezp,
-                                        RealVector& Bxp,
-                                        RealVector& Byp,
-                                        RealVector& Bzp,
-                                        amrex::FArrayBox const * exfab,
-                                        amrex::FArrayBox const * eyfab,
-                                        amrex::FArrayBox const * ezfab,
-                                        amrex::FArrayBox const * bxfab,
-                                        amrex::FArrayBox const * byfab,
-                                        amrex::FArrayBox const * bzfab,
-                                        const int ngE, const int /*e_is_nodal*/,
-                                        const long offset,
-                                        const long np_to_gather,
-                                        int lev,
-                                        int gather_lev)
+PhysicalParticleContainer::PushPX (WarpXParIter& pti,
+                                   amrex::FArrayBox const * exfab,
+                                   amrex::FArrayBox const * eyfab,
+                                   amrex::FArrayBox const * ezfab,
+                                   amrex::FArrayBox const * bxfab,
+                                   amrex::FArrayBox const * byfab,
+                                   amrex::FArrayBox const * bzfab,
+                                   const int ngE, const int /*e_is_nodal*/,
+                                   const long offset,
+                                   const long np_to_push,
+                                   int lev, int gather_lev,
+                                   amrex::Real dt, ScaleFields scaleFields,
+                                   DtType a_dt_type)
 {
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) ||
                                      (gather_lev==(lev  )),
                                      "Gather buffers only work for lev-1");
     // If no particles, do not do anything
     // If do_not_gather = 1 by user, do not do anything
-    if (np_to_gather == 0 || do_not_gather) return;
+    if (np_to_push == 0 || do_not_gather) return;
 
     // Get cell size on gather_lev
     const std::array<Real,3>& dx = WarpX::CellSize(std::max(gather_lev,0));
@@ -1953,6 +1760,7 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti,
     box.grow(ngE);
 
     const auto getPosition = GetParticlePosition(pti, offset);
+          auto setPosition = SetParticlePosition(pti, offset);
 
     const auto getExternalE = GetExternalEField(pti, offset);
     const auto getExternalB = GetExternalBField(pti, offset);
@@ -1966,63 +1774,107 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti,
 
     const Dim3 lo = lbound(box);
 
-    // Depending on l_lower_in_v and WarpX::nox, call
-    // different versions of template function doGatherShapeN
-    if (WarpX::l_lower_order_in_v){
-        if        (WarpX::nox == 1){
-            doGatherShapeN<1,1>(getPosition, getExternalE, getExternalB,
-                                Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
-                                Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
-                                Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
-                                exfab, eyfab, ezfab, bxfab, byfab, bzfab,
-                                np_to_gather, dx,
-                                xyzmin, lo, WarpX::n_rz_azimuthal_modes);
-        } else if (WarpX::nox == 2){
-            doGatherShapeN<2,1>(getPosition, getExternalE, getExternalB,
-                                Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
-                                Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
-                                Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
-                                exfab, eyfab, ezfab, bxfab, byfab, bzfab,
-                                np_to_gather, dx,
-                                xyzmin, lo, WarpX::n_rz_azimuthal_modes);
-        } else if (WarpX::nox == 3){
-            doGatherShapeN<3,1>(getPosition, getExternalE, getExternalB,
-                                Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
-                                Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
-                                Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
-                                exfab, eyfab, ezfab, bxfab, byfab, bzfab,
-                                np_to_gather, dx,
-                                xyzmin, lo, WarpX::n_rz_azimuthal_modes);
-        }
-    } else {
-        if        (WarpX::nox == 1){
-            doGatherShapeN<1,0>(getPosition, getExternalE, getExternalB,
-                                Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
-                                Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
-                                Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
-                                exfab, eyfab, ezfab, bxfab, byfab, bzfab,
-                                np_to_gather, dx,
-                                xyzmin, lo, WarpX::n_rz_azimuthal_modes);
-        } else if (WarpX::nox == 2){
-            doGatherShapeN<2,0>(getPosition, getExternalE, getExternalB,
-                                Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
-                                Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
-                                Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
-                                exfab, eyfab, ezfab, bxfab, byfab, bzfab,
-                                np_to_gather, dx,
-                                xyzmin, lo, WarpX::n_rz_azimuthal_modes);
-        } else if (WarpX::nox == 3){
-            doGatherShapeN<3,0>(getPosition, getExternalE, getExternalB,
-                                Exp.dataPtr() + offset, Eyp.dataPtr() + offset,
-                                Ezp.dataPtr() + offset, Bxp.dataPtr() + offset,
-                                Byp.dataPtr() + offset, Bzp.dataPtr() + offset,
-                                exfab, eyfab, ezfab, bxfab, byfab, bzfab,
-                                np_to_gather, dx,
-                                xyzmin, lo, WarpX::n_rz_azimuthal_modes);
-        }
+    int l_lower_order_in_v = WarpX::l_lower_order_in_v;
+    int nox = WarpX::nox;
+    int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+    amrex::GpuArray<amrex::Real, 3> dx_arr = {dx[0], dx[1], dx[2]};
+    amrex::GpuArray<amrex::Real, 3> xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+    amrex::Array4<const amrex::Real> const& ex_arr = exfab->array();
+    amrex::Array4<const amrex::Real> const& ey_arr = eyfab->array();
+    amrex::Array4<const amrex::Real> const& ez_arr = ezfab->array();
+    amrex::Array4<const amrex::Real> const& bx_arr = bxfab->array();
+    amrex::Array4<const amrex::Real> const& by_arr = byfab->array();
+    amrex::Array4<const amrex::Real> const& bz_arr = bzfab->array();
+
+    amrex::IndexType const ex_type = exfab->box().ixType();
+    amrex::IndexType const ey_type = eyfab->box().ixType();
+    amrex::IndexType const ez_type = ezfab->box().ixType();
+    amrex::IndexType const bx_type = bxfab->box().ixType();
+    amrex::IndexType const by_type = byfab->box().ixType();
+    amrex::IndexType const bz_type = bzfab->box().ixType();
+
+    auto& attribs = pti.GetAttribs();
+    ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr();
+    ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr();
+    ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr();
+
+    auto copyAttribs = CopyParticleAttribs(pti, tmp_particle_data, offset);
+    int do_copy = (WarpX::do_back_transformed_diagnostics &&
+                          do_back_transformed_diagnostics &&
+                   (a_dt_type!=DtType::SecondHalf));
+
+    int* AMREX_RESTRICT ion_lev = nullptr;
+    if (do_field_ionization) {
+        ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
     }
-}
 
+    // Loop over the particles and update their momentum
+    const amrex::Real q = this->charge;
+    const amrex::Real m = this-> mass;
+
+    const auto pusher_algo = WarpX::particle_pusher_algo;
+    const auto do_crr = do_classical_radiation_reaction;
+#ifdef WARPX_QED
+    const auto do_sync = m_do_qed_quantum_sync;
+    amrex::Real t_chi_max = 0.0;
+    if (do_sync) t_chi_max = m_shr_p_qs_engine->get_ref_ctrl().chi_part_min;
+
+    QuantumSynchrotronEvolveOpticalDepth evolve_opt;
+    amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_QSR = nullptr;
+    const bool local_has_quantum_sync = has_quantum_sync();
+    if (local_has_quantum_sync) {
+        evolve_opt = m_shr_p_qs_engine->build_evolve_functor();
+        p_optical_depth_QSR = pti.GetAttribs(particle_comps["optical_depth_QSR"]).dataPtr();
+    }
+#endif
+
+    amrex::ParallelFor( np_to_push, [=] AMREX_GPU_DEVICE (long ip)
+    {
+        amrex::ParticleReal xp, yp, zp;
+        getPosition(ip, xp, yp, zp);
+
+        amrex::ParticleReal Exp = 0._rt, Eyp = 0._rt, Ezp = 0._rt;
+        getExternalE(ip, Exp, Eyp, Ezp);
+
+        amrex::ParticleReal Bxp = 0._rt, Byp = 0._rt, Bzp = 0._rt;
+        getExternalB(ip, Bxp, Byp, Bzp);
+
+        // first gather E and B to the particle positions
+        doGatherShapeN(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                       ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                       ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                       dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes,
+                       nox, l_lower_order_in_v);
+
+        scaleFields(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp);
+
+#ifdef WARPX_WED
+    if (local_has_quantum_sync) {
+        const ParticleReal px = m * ux[ip];
+        const ParticleReal py = m * uy[ip];
+        const ParticleReal pz = m * uz[ip];
+
+        bool has_event_happened = evolve_opt(px, py, pz,
+                                             Exp, Eyp, Ezp,
+                                             Bxp, Byp, Bzp,
+                                             dt, p_optical_depth_QSR[ip]);
+    }
+#endif
+
+        doParticlePush(getPosition, setPosition, copyAttribs, ip,
+                       ux[ip+offset], uy[ip+offset], uz[ip+offset],
+                       Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                       ion_lev ? ion_lev[ip] : 0,
+                       m, q, pusher_algo, do_crr, do_copy,
+#ifdef WARPX_QED
+                       do_sync,
+                       t_chi_max,
+#endif
+                       dt);
+    });
+}
 
 void
 PhysicalParticleContainer::InitIonizationModule ()
@@ -2074,16 +1926,26 @@ PhysicalParticleContainer::InitIonizationModule ()
 }
 
 IonizationFilterFunc
-PhysicalParticleContainer::getIonizationFunc ()
+PhysicalParticleContainer::getIonizationFunc (const WarpXParIter& pti,
+                                              int lev,
+                                              int ngE,
+                                              const amrex::FArrayBox& Ex,
+                                              const amrex::FArrayBox& Ey,
+                                              const amrex::FArrayBox& Ez,
+                                              const amrex::FArrayBox& Bx,
+                                              const amrex::FArrayBox& By,
+                                              const amrex::FArrayBox& Bz)
 {
     WARPX_PROFILE("PPC::getIonizationFunc");
 
-    return IonizationFilterFunc{ionization_energies.dataPtr(),
+    return IonizationFilterFunc(pti, lev, ngE, Ex, Ey, Ez, Bx, By, Bz,
+                                v_galilean,
+                                ionization_energies.dataPtr(),
                                 adk_prefactor.dataPtr(),
                                 adk_exp_prefactor.dataPtr(),
                                 adk_power.dataPtr(),
                                 particle_icomps["ionization_level"],
-                                ion_atomic_number};
+                                ion_atomic_number);
 }
 
 #ifdef WARPX_QED
diff --git a/Source/Particles/Pusher/CopyParticleAttribs.H b/Source/Particles/Pusher/CopyParticleAttribs.H
index 1f8a190ea..b9233d409 100644
--- a/Source/Particles/Pusher/CopyParticleAttribs.H
+++ b/Source/Particles/Pusher/CopyParticleAttribs.H
@@ -17,9 +17,6 @@
 /** \brief Functor that creates copies of the current particle
  *         positions and momenta for later use. This is needed
  *         by the back-transformed diagnostics.
- *
- * \param a_pti iterator to the tile containing the macroparticles
- * \param a_tmp holder for the temporary particle data
 */
 struct CopyParticleAttribs
 {
@@ -39,25 +36,35 @@ struct CopyParticleAttribs
     amrex::ParticleReal* AMREX_RESTRICT uypold = nullptr;
     amrex::ParticleReal* AMREX_RESTRICT uzpold = nullptr;
 
-    CopyParticleAttribs (const WarpXParIter& a_pti, TmpParticles& tmp_particle_data) noexcept
+    /** \brief Construct a new functor
+     *
+     * \param a_pti iterator to the tile containing the macroparticles
+     * \param a_tmp holder for the temporary particle data
+     * \param a_offset offset to apply when reading / writing particle data
+     *        This is needed because when we use field gather buffers we don't
+     *        always start at the particle with index 0.
+     */
+    CopyParticleAttribs (const WarpXParIter& a_pti, TmpParticles& tmp_particle_data,
+                         int a_offset = 0) noexcept
     {
         if (tmp_particle_data.size() == 0) return;
 
         auto& attribs = a_pti.GetAttribs();
-        uxp = attribs[PIdx::ux].dataPtr();
-        uyp = attribs[PIdx::uy].dataPtr();
-        uzp = attribs[PIdx::uz].dataPtr();
+
+        uxp = attribs[PIdx::ux].dataPtr() + a_offset;
+        uyp = attribs[PIdx::uy].dataPtr() + a_offset;
+        uzp = attribs[PIdx::uz].dataPtr() + a_offset;
 
         const auto lev = a_pti.GetLevel();
         const auto index = a_pti.GetPairIndex();
-        xpold  = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr();
-        ypold  = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr();
-        zpold  = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr();
-        uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr();
-        uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr();
-        uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr();
-
-        m_get_position = GetParticlePosition(a_pti);
+        xpold  = tmp_particle_data[lev][index][TmpIdx::xold ].dataPtr() + a_offset;
+        ypold  = tmp_particle_data[lev][index][TmpIdx::yold ].dataPtr() + a_offset;
+        zpold  = tmp_particle_data[lev][index][TmpIdx::zold ].dataPtr() + a_offset;
+        uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr() + a_offset;
+        uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr() + a_offset;
+        uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr() + a_offset;
+
+        m_get_position = GetParticlePosition(a_pti, a_offset);
     }
 
     /** \brief copy the position and momentum of particle i to the
diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H
index cc824ce00..9d4257702 100644
--- a/Source/Particles/RigidInjectedParticleContainer.H
+++ b/Source/Particles/RigidInjectedParticleContainer.H
@@ -68,7 +68,19 @@ public:
                          amrex::Real dt,
                          DtType a_dt_type=DtType::Full) override;
 
-    virtual void PushPX (WarpXParIter& pti, amrex::Real dt, DtType a_dt_type=DtType::Full) override;
+    virtual void PushPX (WarpXParIter& pti,
+                         amrex::FArrayBox const * exfab,
+                         amrex::FArrayBox const * eyfab,
+                         amrex::FArrayBox const * ezfab,
+                         amrex::FArrayBox const * bxfab,
+                         amrex::FArrayBox const * byfab,
+                         amrex::FArrayBox const * bzfab,
+                         const int ngE, const int /*e_is_nodal*/,
+                         const long offset,
+                         const long np_to_push,
+                         int lev, int gather_lev,
+                         amrex::Real dt, ScaleFields scaleFields,
+                         DtType a_dt_type=DtType::Full) override;
 
     virtual void PushP (int lev, amrex::Real dt,
                         const amrex::MultiFab& Ex,
diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp
index ee0890fd7..6d5565f2e 100644
--- a/Source/Particles/RigidInjectedParticleContainer.cpp
+++ b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -20,6 +20,8 @@
 #include "Pusher/UpdateMomentumBorisWithRadiationReaction.H"
 #include "Pusher/UpdateMomentumHigueraCary.H"
 #include "Pusher/GetAndSetPosition.H"
+#include "Gather/ScaleFields.H"
+#include "Gather/FieldGather.H"
 
 #include <limits>
 #include <sstream>
@@ -224,10 +226,20 @@ RigidInjectedParticleContainer::BoostandRemapParticles()
 }
 
 void
-RigidInjectedParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_type)
+RigidInjectedParticleContainer::PushPX (WarpXParIter& pti,
+                                        amrex::FArrayBox const * exfab,
+                                        amrex::FArrayBox const * eyfab,
+                                        amrex::FArrayBox const * ezfab,
+                                        amrex::FArrayBox const * bxfab,
+                                        amrex::FArrayBox const * byfab,
+                                        amrex::FArrayBox const * bzfab,
+                                        const int ngE, const int e_is_nodal,
+                                        const long offset,
+                                        const long np_to_push,
+                                        int lev, int gather_lev,
+                                        amrex::Real dt, ScaleFields /*scaleFields*/,
+                                        DtType a_dt_type)
 {
-
-    // This wraps the momentum and position advance so that inheritors can modify the call.
     auto& attribs = pti.GetAttribs();
     auto& uxp = attribs[PIdx::ux];
     auto& uyp = attribs[PIdx::uy];
@@ -243,12 +255,6 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_
     ParticleReal* const AMREX_RESTRICT ux = uxp.dataPtr();
     ParticleReal* const AMREX_RESTRICT uy = uyp.dataPtr();
     ParticleReal* const AMREX_RESTRICT uz = uzp.dataPtr();
-    ParticleReal* const AMREX_RESTRICT Exp = attribs[PIdx::Ex].dataPtr();
-    ParticleReal* const AMREX_RESTRICT Eyp = attribs[PIdx::Ey].dataPtr();
-    ParticleReal* const AMREX_RESTRICT Ezp = attribs[PIdx::Ez].dataPtr();
-    ParticleReal* const AMREX_RESTRICT Bxp = attribs[PIdx::Bx].dataPtr();
-    ParticleReal* const AMREX_RESTRICT Byp = attribs[PIdx::By].dataPtr();
-    ParticleReal* const AMREX_RESTRICT Bzp = attribs[PIdx::Bz].dataPtr();
 
     if (!done_injecting_lev)
     {
@@ -282,32 +288,15 @@ RigidInjectedParticleContainer::PushPX (WarpXParIter& pti, Real dt, DtType a_dt_
                                 uyp_save_ptr[i] = uy[i];
                                 uzp_save_ptr[i] = uz[i];
                             });
-
-        // Scale the fields of particles about to cross the injection plane.
-        // This only approximates what should be happening. The particles
-        // should by advanced a fraction of a time step instead.
-        // Scaling the fields is much easier and may be good enough.
-        const Real v_boost = WarpX::beta_boost*PhysConst::c;
-        const Real z_plane_previous = zinject_plane_lev_previous;
-        const Real vz_ave_boosted = vzbeam_ave_boosted;
-        amrex::ParallelFor( pti.numParticles(),
-            [=] AMREX_GPU_DEVICE (long i) {
-                                ParticleReal xp, yp, zp;
-                                GetPosition(i, xp, yp, zp);
-                                const Real dtscale = dt - (z_plane_previous - zp)/(vz_ave_boosted + v_boost);
-                                if (0. < dtscale && dtscale < dt) {
-                                    Exp[i] *= dtscale;
-                                    Eyp[i] *= dtscale;
-                                    Ezp[i] *= dtscale;
-                                    Bxp[i] *= dtscale;
-                                    Byp[i] *= dtscale;
-                                    Bzp[i] *= dtscale;
-                                }
-                            }
-            );
     }
 
-    PhysicalParticleContainer::PushPX(pti, dt, a_dt_type);
+    const bool do_scale = not done_injecting_lev;
+    const Real v_boost = WarpX::beta_boost*PhysConst::c;
+    PhysicalParticleContainer::PushPX(pti, exfab, eyfab, ezfab, bxfab, byfab, bzfab,
+                                      ngE, e_is_nodal, offset, np_to_push, lev, gather_lev, dt,
+                                      ScaleFields(do_scale, dt, zinject_plane_lev_previous,
+                                                  vzbeam_ave_boosted, v_boost),
+                                      a_dt_type);
 
     if (!done_injecting_lev) {
 
@@ -395,23 +384,16 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
 
     if (do_not_push) return;
 
+    const std::array<Real,3>& dx = WarpX::CellSize(std::max(lev,0));
+
 #ifdef _OPENMP
 #pragma omp parallel
 #endif
     {
         for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
         {
-            auto& attribs = pti.GetAttribs();
-
-            auto& uxp = attribs[PIdx::ux];
-            auto& uyp = attribs[PIdx::uy];
-            auto& uzp = attribs[PIdx::uz];
-            auto& Exp = attribs[PIdx::Ex];
-            auto& Eyp = attribs[PIdx::Ey];
-            auto& Ezp = attribs[PIdx::Ez];
-            auto& Bxp = attribs[PIdx::Bx];
-            auto& Byp = attribs[PIdx::By];
-            auto& Bzp = attribs[PIdx::Bz];
+            amrex::Box box = pti.tilebox();
+            box.grow(Ex.nGrow());
 
             const long np = pti.numParticles();
 
@@ -423,76 +405,108 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
             const FArrayBox& byfab = By[pti];
             const FArrayBox& bzfab = Bz[pti];
 
-            int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
-            FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
-                        Ex.nGrow(), e_is_nodal,
-                        0, np, lev, lev);
+            const auto getPosition = GetParticlePosition(pti);
+                  auto setPosition = SetParticlePosition(pti);
+
+            const auto getExternalE = GetExternalEField(pti);
+            const auto getExternalB = GetExternalBField(pti);
+
+            const auto& xyzmin = WarpX::GetInstance().LowerCornerWithGalilean(box,v_galilean,lev);
+
+            const Dim3 lo = lbound(box);
+
+            int l_lower_order_in_v = WarpX::l_lower_order_in_v;
+            int nox = WarpX::nox;
+            int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+
+            amrex::GpuArray<amrex::Real, 3> dx_arr = {dx[0], dx[1], dx[2]};
+            amrex::GpuArray<amrex::Real, 3> xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]};
+
+            amrex::Array4<const amrex::Real> const& ex_arr = exfab.array();
+            amrex::Array4<const amrex::Real> const& ey_arr = eyfab.array();
+            amrex::Array4<const amrex::Real> const& ez_arr = ezfab.array();
+            amrex::Array4<const amrex::Real> const& bx_arr = bxfab.array();
+            amrex::Array4<const amrex::Real> const& by_arr = byfab.array();
+            amrex::Array4<const amrex::Real> const& bz_arr = bzfab.array();
+
+            amrex::IndexType const ex_type = exfab.box().ixType();
+            amrex::IndexType const ey_type = eyfab.box().ixType();
+            amrex::IndexType const ez_type = ezfab.box().ixType();
+            amrex::IndexType const bx_type = bxfab.box().ixType();
+            amrex::IndexType const by_type = byfab.box().ixType();
+            amrex::IndexType const bz_type = bzfab.box().ixType();
+
+            auto& attribs = pti.GetAttribs();
+            auto& uxp = attribs[PIdx::ux];
+            auto& uyp = attribs[PIdx::uy];
+            auto& uzp = attribs[PIdx::uz];
+            amrex::ParticleReal* const AMREX_RESTRICT uxpp = attribs[PIdx::ux].dataPtr();
+            amrex::ParticleReal* const AMREX_RESTRICT uypp = attribs[PIdx::uy].dataPtr();
+            amrex::ParticleReal* const AMREX_RESTRICT uzpp = attribs[PIdx::uz].dataPtr();
+
+            int* AMREX_RESTRICT ion_lev = nullptr;
+            if (do_field_ionization) {
+                ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
+            }
 
             // Save the position and momenta, making copies
             auto uxp_save = uxp;
             auto uyp_save = uyp;
             auto uzp_save = uzp;
 
-            // This wraps the momentum advance so that inheritors can modify the call.
-            // Extract pointers to the different particle quantities
-            const auto GetPosition = GetParticlePosition(pti);
-            ParticleReal* const AMREX_RESTRICT uxpp = uxp.dataPtr();
-            ParticleReal* const AMREX_RESTRICT uypp = uyp.dataPtr();
-            ParticleReal* const AMREX_RESTRICT uzpp = uzp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Expp = Exp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Eypp = Eyp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Ezpp = Ezp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Bxpp = Bxp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Bypp = Byp.dataPtr();
-            const ParticleReal* const AMREX_RESTRICT Bzpp = Bzp.dataPtr();
-
             // Loop over the particles and update their momentum
-            const Real q = this->charge;
-            const Real m = this->mass;
-
-            //Assumes that all consistency checks have been done at initialization
-            if(do_classical_radiation_reaction){
-                amrex::ParallelFor(
-                    pti.numParticles(),
-                    [=] AMREX_GPU_DEVICE (long i) {
-                        UpdateMomentumBorisWithRadiationReaction(
-                        uxpp[i], uypp[i], uzpp[i],
-                        Expp[i], Eypp[i], Ezpp[i],
-                        Bxpp[i], Bypp[i], Bzpp[i],
-                        q, m, dt);
-                    }
-                );
-            } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
-                amrex::ParallelFor( pti.numParticles(),
-                    [=] AMREX_GPU_DEVICE (long i) {
-                        UpdateMomentumBoris(
-                            uxpp[i], uypp[i], uzpp[i],
-                            Expp[i], Eypp[i], Ezpp[i],
-                            Bxpp[i], Bypp[i], Bzpp[i],
-                            q, m, dt);
-                    }
-                );
-            } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) {
-                amrex::ParallelFor( pti.numParticles(),
-                    [=] AMREX_GPU_DEVICE (long i) {
-                        UpdateMomentumVay(
-                            uxpp[i], uypp[i], uzpp[i],
-                            Expp[i], Eypp[i], Ezpp[i],
-                            Bxpp[i], Bypp[i], Bzpp[i],
-                            q, m, dt);
-                    }
-                );
-            } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) {
-                amrex::ParallelFor( pti.numParticles(),
-                    [=] AMREX_GPU_DEVICE (long i) {
-                        UpdateMomentumHigueraCary( uxpp[i], uypp[i], uzpp[i],
-                              Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
-                    }
-                );
-            } else {
-              amrex::Abort("Unknown particle pusher");
-            }
+            const amrex::Real q = this->charge;
+            const amrex::Real m = this-> mass;
+
+            const auto pusher_algo = WarpX::particle_pusher_algo;
+            const auto do_crr = do_classical_radiation_reaction;
+
+            amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip)
+            {
+                amrex::ParticleReal xp, yp, zp;
+                getPosition(ip, xp, yp, zp);
+
+                amrex::ParticleReal Exp = 0._rt, Eyp = 0._rt, Ezp = 0._rt;
+                getExternalE(ip, Exp, Eyp, Ezp);
+
+                amrex::ParticleReal Bxp = 0._rt, Byp = 0._rt, Bzp = 0._rt;
+                getExternalB(ip, Bxp, Byp, Bzp);
+
+                // first gather E and B to the particle positions
+                doGatherShapeN(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
+                               ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr,
+                               ex_type, ey_type, ez_type, bx_type, by_type, bz_type,
+                               dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes,
+                               nox, l_lower_order_in_v);
+
+                if (do_crr) {
+                    amrex::Real qp = q;
+                    if (ion_lev) { qp *= ion_lev[ip]; }
+                    UpdateMomentumBorisWithRadiationReaction(uxpp[ip], uypp[ip], uzpp[ip],
+                                                             Exp, Eyp, Ezp, Bxp,
+                                                             Byp, Bzp, qp, m, dt);
+                } else if (pusher_algo == ParticlePusherAlgo::Boris) {
+                    amrex::Real qp = q;
+                    if (ion_lev) { qp *= ion_lev[ip]; }
+                    UpdateMomentumBoris( uxpp[ip], uypp[ip], uzpp[ip],
+                                         Exp, Eyp, Ezp, Bxp,
+                                         Byp, Bzp, qp, m, dt);
+                } else if (pusher_algo == ParticlePusherAlgo::Vay) {
+                    amrex::Real qp = q;
+                    if (ion_lev){ qp *= ion_lev[ip]; }
+                    UpdateMomentumVay( uxpp[ip], uypp[ip], uzpp[ip],
+                                       Exp, Eyp, Ezp, Bxp,
+                                       Byp, Bzp, qp, m, dt);
+                } else if (pusher_algo == ParticlePusherAlgo::HigueraCary) {
+                    amrex::Real qp = q;
+                    if (ion_lev){ qp *= ion_lev[ip]; }
+                    UpdateMomentumHigueraCary( uxpp[ip], uypp[ip], uzpp[ip],
+                                               Exp, Eyp, Ezp, Bxp,
+                                               Byp, Bzp, qp, m, dt);
+                } else {
+                    amrex::Abort("Unknown particle pusher");
+                }
+            });
 
             // Undo the push for particles not injected yet.
             // It is assumed that PushP will only be called on the first and last steps
@@ -501,17 +515,16 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
             const ParticleReal* const AMREX_RESTRICT uy_save = uyp_save.dataPtr();
             const ParticleReal* const AMREX_RESTRICT uz_save = uzp_save.dataPtr();
             const ParticleReal zz = zinject_plane_levels[lev];
-            amrex::ParallelFor( pti.numParticles(),
-                                [=] AMREX_GPU_DEVICE (long i) {
-                                    ParticleReal xp, yp, zp;
-                                    GetPosition(i, xp, yp, zp);
-                                    if (zp <= zz) {
-                                        uxpp[i] = ux_save[i];
-                                        uypp[i] = uy_save[i];
-                                        uzpp[i] = uz_save[i];
-                                    }
-                                }
-                );
+            amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i)
+            {
+                ParticleReal xp, yp, zp;
+                getPosition(i, xp, yp, zp);
+                if (zp <= zz) {
+                    uxpp[i] = ux_save[i];
+                    uypp[i] = uy_save[i];
+                    uzpp[i] = uz_save[i];
+                }
+            });
         }
     }
 }
diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H
index 8bc4491d7..33769c61c 100644
--- a/Source/Particles/WarpXParticleContainer.H
+++ b/Source/Particles/WarpXParticleContainer.H
@@ -32,7 +32,7 @@ struct PIdx
 {
     enum { // Particle Attributes stored in amrex::ParticleContainer's struct of array
         w = 0,  // weight
-        ux, uy, uz, Ex, Ey, Ez, Bx, By, Bz,
+        ux, uy, uz,
 #ifdef WARPX_DIM_RZ
         theta, // RZ needs all three position components
 #endif
@@ -65,14 +65,8 @@ namespace ParticleStringNames
         {"ux",    PIdx::ux   },
         {"uy",    PIdx::uy   },
         {"uz",    PIdx::uz   },
-        {"Ex",    PIdx::Ex   },
-        {"Ey",    PIdx::Ey   },
-        {"Ez",    PIdx::Ez   },
-        {"Bx",    PIdx::Bx   },
-        {"By",    PIdx::By   },
-        {"Bz",    PIdx::Bz   }
 #ifdef WARPX_DIM_RZ
-        ,{"theta", PIdx::theta}
+        {"theta", PIdx::theta}
 #endif
     };
 }
@@ -85,6 +79,8 @@ public:
 
     WarpXParIter (ContainerType& pc, int level);
 
+    WarpXParIter (ContainerType& pc, int level, amrex::MFItInfo& info);
+
     const std::array<RealVector, PIdx::nattribs>& GetAttribs () const {
         return GetStructOfArrays().GetRealData();
     }
@@ -128,9 +124,7 @@ class MultiParticleContainer;
  * push.
  *
  * Note: many functions are pure virtual (meaning they MUST be defined in
- * derived classes, e.g., Evolve) or empty function (meaning they
- * do not do anything, e.g., FieldGather, meant to be overriden by derived
- * function) or actual functions (e.g. CurrentDeposition).
+ * derived classes, e.g., Evolve) or actual functions (e.g. CurrentDeposition).
  */
 class WarpXParticleContainer
     : public amrex::ParticleContainer<0,0,PIdx::nattribs>
@@ -155,14 +149,6 @@ public:
 
     virtual void InitData () = 0;
 
-    virtual void FieldGatherES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E,
-                                const amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<int> > > >& masks) {}
-
-    virtual void FieldGather (int lev,
-                              const amrex::MultiFab& Ex, const amrex::MultiFab& Ey,
-                              const amrex::MultiFab& Ez, const amrex::MultiFab& Bx,
-                              const amrex::MultiFab& By, const amrex::MultiFab& Bz) {}
-
     /**
      * Evolve is the central WarpXParticleContainer function that advances
      * particles for a time dt (typically one timestep). It is a pure virtual
@@ -342,6 +328,8 @@ public:
     template<PhysicalSpecies PhysSpec>
     bool AmIA () const noexcept {return (physical_species == PhysSpec);}
 
+    amrex::Array<amrex::Real,3> get_v_galilean () {return v_galilean;}
+
 protected:
     amrex::Array<amrex::Real,3> v_galilean = {{0}};
     std::map<std::string, int> particle_comps;
diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp
index 27d042f74..95e656a4f 100644
--- a/Source/Particles/WarpXParticleContainer.cpp
+++ b/Source/Particles/WarpXParticleContainer.cpp
@@ -26,7 +26,14 @@
 using namespace amrex;
 
 WarpXParIter::WarpXParIter (ContainerType& pc, int level)
-    : ParIter(pc, level, MFItInfo().SetDynamic(WarpX::do_dynamic_scheduling))
+    : amrex::ParIter<0,0,PIdx::nattribs>(pc, level,
+             MFItInfo().SetDynamic(WarpX::do_dynamic_scheduling))
+{
+}
+
+WarpXParIter::WarpXParIter (ContainerType& pc, int level, MFItInfo& info)
+    : amrex::ParIter<0,0,PIdx::nattribs>(pc, level,
+                   info.SetDynamic(WarpX::do_dynamic_scheduling))
 {
 }
 
@@ -34,9 +41,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies)
     : ParticleContainer<0,0,PIdx::nattribs>(amr_core->GetParGDB())
     , species_id(ispecies)
 {
-    for (unsigned int i = PIdx::Ex; i <= PIdx::Bz; ++i) {
-        communicate_real_comp[i] = false; // Don't need to communicate E and B.
-    }
     SetParticleSize();
     ReadParameters();
 
@@ -45,12 +49,6 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies)
     particle_comps["ux"] = PIdx::ux;
     particle_comps["uy"] = PIdx::uy;
     particle_comps["uz"] = PIdx::uz;
-    particle_comps["Ex"] = PIdx::Ex;
-    particle_comps["Ey"] = PIdx::Ey;
-    particle_comps["Ez"] = PIdx::Ez;
-    particle_comps["Bx"] = PIdx::Bx;
-    particle_comps["By"] = PIdx::By;
-    particle_comps["Bz"] = PIdx::Bz;
 #ifdef WARPX_DIM_RZ
     particle_comps["theta"] = PIdx::theta;
 #endif
diff --git a/Source/WarpX.H b/Source/WarpX.H
index fe00115fd..e080767b7 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -292,10 +292,10 @@ public:
     void EvolveB (int lev, PatchType patch_type, amrex::Real dt);
     void EvolveE (int lev, PatchType patch_type, amrex::Real dt);
     void EvolveF (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type);
+
     void MacroscopicEvolveE (         amrex::Real dt);
     void MacroscopicEvolveE (int lev, amrex::Real dt);
     void MacroscopicEvolveE (int lev, PatchType patch_type, amrex::Real dt);
-    void FieldGather ();
 
     /** \brief apply QED correction on electric field
      * \param dt vector of time steps (for all levels)
@@ -350,6 +350,22 @@ public:
 
     PML* GetPML (int lev);
 
+    /** Run the ionization module on all species */
+    void doFieldIonization ();
+    /** Run the ionization module on all species at level lev
+     * \param lev level
+     */
+    void doFieldIonization (int lev);
+
+#ifdef WARPX_QED
+    /** Run the QED module on all species */
+    void doQEDEvents ();
+    /** Run the QED module on all species at level lev
+     * \param lev level
+     */
+    void doQEDEvents (int lev);
+#endif
+
     void PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type=DtType::Full);
     void PushParticlesandDepose (         amrex::Real cur_time);
 
@@ -405,6 +421,11 @@ public:
                                                   std::array<amrex::Real,3> galilean_shift, int lev);
     static std::array<amrex::Real,3> UpperCorner (const amrex::Box& bx, int lev);
 
+    /*
+      /brief This computes the lower of the problem domain, taking into account any shift when using the Galilean algorithm.
+     */
+    std::array<amrex::Real,3> LowerCornerWithGalilean (const amrex::Box& bx, const amrex::Array<amrex::Real,3>& v_galilean, int lev);
+
     static amrex::IntVect RefRatio (int lev);
 
     static const amrex::iMultiFab* CurrentBufferMasks (int lev);
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 476474b9d..8139fc92a 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -1102,6 +1102,15 @@ WarpX::UpperCorner(const Box& bx, int lev)
 #endif
 }
 
+std::array<Real,3>
+WarpX::LowerCornerWithGalilean (const Box& bx, const amrex::Array<amrex::Real,3>& v_galilean, int lev)
+{
+    amrex::Real cur_time = gett_new(lev);
+    amrex::Real time_shift = (cur_time - time_of_last_gal_shift);
+    amrex::Array<amrex::Real,3> galilean_shift = { v_galilean[0]*time_shift, v_galilean[1]*time_shift, v_galilean[2]*time_shift };
+    return WarpX::LowerCorner(bx, galilean_shift, lev);
+}
+
 IntVect
 WarpX::RefRatio (int lev)
 {
@@ -1361,13 +1370,3 @@ WarpX::PicsarVersion ()
     return std::string("Unknown");
 #endif
 }
-
-void
-WarpX::FieldGather ()
-{
-    for (int lev = 0; lev <= finest_level; ++lev) {
-        mypc->FieldGather(lev,
-                          *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2],
-                          *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]);
-    }
-}
author	Andrew Myers <atmyers@lbl.gov>	2020-06-22 09:28:16 -0700
committer	GitHub <noreply@github.com>	2020-06-22 09:28:16 -0700
commit	f4978e1001494e2b148128380fb37f3c2450209f (patch)
tree	237a671ca194338411f16f52797e70114b220384
parent	97755d1c2e04e5b8c3295182eaff472c73cf8a53 (diff)
download	WarpX-f4978e1001494e2b148128380fb37f3c2450209f.tar.gz WarpX-f4978e1001494e2b148128380fb37f3c2450209f.tar.zst WarpX-f4978e1001494e2b148128380fb37f3c2450209f.zip