diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 77d1ac542cf..f254f9720b9 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -4271,14 +4271,31 @@ In-situ capabilities can be used by turning on Sensei or Ascent (provided they a
     List of particle quantities to write to output.
     Choices are ``x``, ``y``, ``z`` for the particle positions (3D, RZ, RSPHERE), ``x`` and ``z`` in 2D, ``z`` in 1D, ``x`` and ``y`` for RCYLINDER,
     ``w`` for the particle weight and ``ux``, ``uy``, ``uz`` for the particle momenta.
-    When using the lab-frame electrostatic solver, ``phi`` (electrostatic potential, on the macroparticles) is also available.
-    By default, all particle quantities (except ``phi``) are written.
+    When writing to the OpenPMD format, the fields can also be obtained, ``Ex``, ``Ey``, ``Ez``, ``Bx``, ``By``, ``Bz``.
+    Also, when writing to the OpenPMD format and when using the lab-frame electrostatic solver, ``phi`` (electrostatic potential, on the macroparticles) is also available.
+    By default, positions, momenta, and weights are written out.
     If ``<diag_name>.<species_name>.variables = none``, no particle data are written.
 
+    .. note::
+        The electromagnetic fields gathered on the particles are computed in the same way as during the simulation
+        (i.e. interpolating ``Efield_aux`` / ``Bfield_aux`` with the simulation's particle shape factor),
+        and do not include any externally applied fields.
+
+.. pp:param:: <diag_name>.<species_name>.additional_variables
+    :type: list of `strings` separated by spaces
+    :optional:
+
+    List of additional particle quantities to write to output, when using the OpenPMD format.
+    This allows specifying the additional particle quantities beyond the standard position, momentum, and weight.
+    The options are the fields, ``Ex``, ``Ey``, ``Ez``, ``Bx``, ``By``, ``Bz``,
+    and when using the lab-frame electrostatic solver, the electrostatic potential ``phi``.
+    As for ``variables``, the gathered electromagnetic fields do not include any externally applied fields.
+
 .. pp:param:: <diag_name>.<species_name>.random_fraction
     :type: `float`
     :optional:
 
+
     If provided ``<diag_name>.<species_name>.random_fraction = a``, only `a` fraction of the particle data of this species will be dumped randomly in diag ``<diag_name>``, i.e. if `rand() < a`, this particle will be dumped, where `rand()` denotes a random number generator.
     The value `a` provided should be between 0 and 1.
 
diff --git a/Examples/Tests/CMakeLists.txt b/Examples/Tests/CMakeLists.txt
index bf5a4ddd897..aabc394bd45 100644
--- a/Examples/Tests/CMakeLists.txt
+++ b/Examples/Tests/CMakeLists.txt
@@ -19,6 +19,7 @@ add_subdirectory(embedded_boundary_em_particle_absorption)
 add_subdirectory(embedded_boundary_python_api)
 add_subdirectory(embedded_boundary_rotated_cube)
 add_subdirectory(embedded_circle)
+add_subdirectory(EM_fields_on_particles)
 add_subdirectory(energy_conserving_thermal_plasma)
 add_subdirectory(field_probe)
 add_subdirectory(flux_injection)
diff --git a/Examples/Tests/EM_fields_on_particles/CMakeLists.txt b/Examples/Tests/EM_fields_on_particles/CMakeLists.txt
new file mode 100644
index 00000000000..4ef51e5894e
--- /dev/null
+++ b/Examples/Tests/EM_fields_on_particles/CMakeLists.txt
@@ -0,0 +1,11 @@
+# add tests
+
+add_warpx_test(
+    test_2d_particle_EM_diagnostics # name
+    2 # dims
+    2 # nprocs
+    inputs_test_2d_particle_EM_diagnostics # inputs
+    "analysis_2d.py" # analysis
+    "analysis_default_regression.py --path diags/diag_checksum/" # checksum
+    OFF # dependency
+)
diff --git a/Examples/Tests/EM_fields_on_particles/analysis_2d.py b/Examples/Tests/EM_fields_on_particles/analysis_2d.py
new file mode 100755
index 00000000000..9c1c5f8ee6e
--- /dev/null
+++ b/Examples/Tests/EM_fields_on_particles/analysis_2d.py
@@ -0,0 +1,165 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+from scipy.interpolate import interp1d
+from scipy.optimize import minimize
+from scipy.signal import hilbert
+
+c = 3e8
+lambda_laser = 800e-9
+T_laser = lambda_laser / c
+
+E_max = 3.26e10
+T_peak = 10e-15
+tau = 5e-15
+delta_t = 1.6e-6 / c  # time for the laser to reach the particle
+
+ts_particle = OpenPMDTimeSeries("diags/diag1/")
+
+ex_t, ey_t, ez_t, bx_t, by_t, bz_t = ts_particle.iterate(
+    ts_particle.get_particle, ["ex", "ey", "ez", "bx", "by", "bz"], species="electron"
+)
+
+ex_t = np.squeeze(ex_t)
+ey_t = np.squeeze(ey_t)
+ez_t = np.squeeze(ez_t)
+bx_t = np.squeeze(bx_t)
+by_t = np.squeeze(by_t)
+bz_t = np.squeeze(bz_t)
+
+n_diags = 401
+
+DT = 6.324524234e-17  # @ CFL = 0.99
+iterations = np.linspace(0, n_diags - 1, n_diags)
+T = (iterations * DT) - DT / 2
+
+T_peak_part = T_peak + delta_t
+Ey_max = E_max * 1 / np.sqrt(5)  # polarisation vector (0 1 2)
+Ez_max = E_max * np.sqrt(1 - 1 / 5)
+By_max = Ez_max / c
+Bz_max = Ey_max / c
+
+
+def get_laser_th(E_0, T_p, T, tau, T_laser, phi=0):
+    alpha = np.exp(-((T - T_peak_part) ** 2) / (tau**2))
+    return E_0 * alpha * np.cos(2 * np.pi * (T - T_p) / T_laser + phi)
+
+
+Ey_th = get_laser_th(Ey_max, T_peak_part, T, tau, T_laser)
+Ez_th = get_laser_th(Ez_max, T_peak_part, T, tau, T_laser)
+By_th = get_laser_th(By_max, T_peak_part, T, tau, T_laser, phi=np.pi)
+Bz_th = get_laser_th(Bz_max, T_peak_part, T, tau, T_laser)
+
+Ey_th = np.where(T >= delta_t, Ey_th, 0)
+Ez_th = np.where(T >= delta_t, Ez_th, 0)
+By_th = np.where(T >= delta_t, By_th, 0)
+Bz_th = np.where(T >= delta_t, Bz_th, 0)
+
+
+# due to injection method the simulated fields have a slight dephasing
+# wrt the theoretical field
+# here we compute the dephasing and correct it on the theoretical field
+def align_cost(shift, sig1, sig2, time):
+    # Create an interpolation function for the signal to be shifted
+    interp_func = interp1d(time, sig1, kind="linear", fill_value="extrapolate")
+    shifted_sig1 = interp_func(time + shift)
+    # MSE
+    cost = np.mean((shifted_sig1 - sig2) ** 2)
+    return cost
+
+
+init_shift = 0.0
+res = minimize(
+    align_cost,
+    init_shift,
+    args=(Ey_th, ey_t, T),
+    method="Nelder-Mead",
+    options={"xatol": 1e-20, "fatol": 1e-20, "maxiter": 10000},
+)
+opt_shift = res.x[0]
+
+interp_func_ey = interp1d(T, Ey_th, kind="cubic", fill_value="extrapolate")
+interp_func_ez = interp1d(T, Ez_th, kind="cubic", fill_value="extrapolate")
+interp_func_by = interp1d(T, By_th, kind="cubic", fill_value="extrapolate")
+interp_func_bz = interp1d(T, Bz_th, kind="cubic", fill_value="extrapolate")
+
+Ey_aligned = interp_func_ey(T + opt_shift)
+Ez_aligned = interp_func_ez(T + opt_shift)
+By_aligned = interp_func_by(T + opt_shift)
+Bz_aligned = interp_func_bz(T + opt_shift)
+
+# test cases
+
+# E.1
+
+assert np.allclose(np.zeros(ex_t.shape), ex_t, rtol=0, atol=1e-2)
+
+# E.2
+
+ey_p = ey_t[T >= delta_t]
+Ey_p = Ey_aligned[T >= delta_t]
+ey_p += 1e11
+Ey_p += 1e11
+
+# np.allclose will almost always be false when values are very close to zero
+# unless we give a high absolute tolerance
+# Here I offset the values by 1e11 to avoid this problem
+
+ez_p = ez_t[T >= delta_t]
+Ez_p = Ez_aligned[T >= delta_t]
+ez_p += 1e11
+Ez_p += 1e11
+
+assert np.allclose(ey_p, Ey_p, rtol=8e-4, atol=0)
+assert np.allclose(ez_p, Ez_p, rtol=2e-3, atol=0)
+
+# E.3
+
+ey_m = ey_t[T < delta_t - DT * 4]
+ez_m = ez_t[T < delta_t - DT * 4]
+assert np.allclose(ey_m, np.zeros(ey_m.shape), atol=1e-10, rtol=0)
+assert np.allclose(ez_m, np.zeros(ey_m.shape), atol=1e-10, rtol=0)
+
+# E.4
+
+env_ez_t = np.abs(hilbert(ez_t))
+env_ey_t = np.abs(hilbert(ey_t))
+assert np.isclose(np.max(env_ez_t) / np.max(env_ey_t), 2, rtol=2e-3, atol=0)
+assert np.allclose(2 * env_ey_t[85:390], env_ez_t[85:390], rtol=8e-3, atol=0)
+
+# B.1
+
+assert np.allclose(np.zeros(bx_t.shape), bx_t, rtol=0, atol=1e-3)
+
+# B.2
+
+by_p = by_t[T >= delta_t]
+By_p = By_aligned[T >= delta_t]
+by_p += 1e3
+By_p += 1e3
+
+# same as for E.2
+# np.allclose will almost always be false when values are very close to zero
+
+bz_p = bz_t[T >= delta_t]
+Bz_p = Bz_aligned[T >= delta_t]
+bz_p += 1e3
+Bz_p += 1e3
+
+assert np.allclose(by_p, By_p, rtol=8e-4, atol=0)
+assert np.allclose(bz_p, Bz_p, rtol=4e-4, atol=0)
+
+# B.3
+
+by_m = by_t[T < delta_t - DT * 4]
+bz_m = bz_t[T < delta_t - DT * 4]
+assert np.allclose(by_m, np.zeros(by_m.shape), atol=1e-10, rtol=0)
+assert np.allclose(bz_m, np.zeros(by_m.shape), atol=1e-10, rtol=0)
+
+# B.4
+
+env_bz_t = np.abs(hilbert(bz_t))
+env_by_t = np.abs(hilbert(by_t))
+assert np.isclose(np.max(env_by_t) / np.max(env_bz_t), 2, rtol=2e-3, atol=0)
+assert np.allclose(0.5 * env_by_t[85:390], env_bz_t[85:390], rtol=8e-3, atol=0)
diff --git a/Examples/Tests/EM_fields_on_particles/analysis_default_regression.py b/Examples/Tests/EM_fields_on_particles/analysis_default_regression.py
new file mode 120000
index 00000000000..d8ce3fca419
--- /dev/null
+++ b/Examples/Tests/EM_fields_on_particles/analysis_default_regression.py
@@ -0,0 +1 @@
+../../analysis_default_regression.py
\ No newline at end of file
diff --git a/Examples/Tests/EM_fields_on_particles/inputs_test_2d_particle_EM_diagnostics b/Examples/Tests/EM_fields_on_particles/inputs_test_2d_particle_EM_diagnostics
new file mode 100644
index 00000000000..a13df22d624
--- /dev/null
+++ b/Examples/Tests/EM_fields_on_particles/inputs_test_2d_particle_EM_diagnostics
@@ -0,0 +1,84 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+
+# Parameters
+max_step = 400
+amr.n_cell = 256 256
+warpx.numprocs = 2 1
+amr.max_level = 0
+geometry.dims = 2
+geometry.prob_lo = 0 0   # physical domain
+geometry.prob_hi = 5e-6 2.5e-5
+
+# BC
+
+boundary.field_lo = pml pml
+boundary.field_hi = pml pml
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.serialize_initial_conditions = 1
+warpx.verbose = 1
+warpx.cfl = 0.99
+warpx.use_filter = 0
+
+# Order of particle shape factors
+algo.particle_shape = 2
+
+#################################
+########### "PLASMA" ############
+#################################
+
+particles.species_names = electron
+
+electron.species_type = electron
+electron.injection_style = "SingleParticle"
+electron.single_particle_pos = 2.5e-6 0 1.25e-5
+electron.single_particle_u = 0 0 0
+electron.single_particle_weight = 0
+
+electron.do_not_push = 1
+electron.do_not_deposit = 1
+
+#################################
+######### LASER #################
+#################################
+
+lasers.names = las1
+
+las1.position = 0.9e-6 326 1.25e-5
+las1.polarization = 0 1 2
+las1.direction = 1 0 0
+las1.e_max = 3.26e10
+las1.wavelength = 800e-9
+las1.profile = "Gaussian"
+las1.profile_t_peak = 10e-15
+las1.profile_duration = 5e-15
+las1.profile_waist = 5e-6
+las1.profile_focal_distance = 1.6e-6
+
+#################################
+############# DIAGS #############
+#################################
+
+diagnostics.diags_names = diag1 diag_checksum
+
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.format = openpmd
+diag1.openpmd_backend = h5
+diag1.fields_to_plot = none
+diag1.write_species = 1
+diag1.species = electron
+diag1.electron.variables = Ex Ey Ez Bx By Bz # will emit a warning because no position is asked
+
+diag_checksum.intervals = 400:400
+diag_checksum.diag_type = Full
+diag_checksum.format = openpmd
+diag_checksum.openpmd_backend = h5
+diag_checksum.fields_to_plot = Ex Ey Ez Bx By Bz
+diag_checksum.write_species = 1
+diag_checksum.species = electron
+diag_checksum.electron.variables = x z ux uy uz w Ex Ey Ez Bx By Bz
diff --git a/Examples/Tests/langmuir/analysis_3d.py b/Examples/Tests/langmuir/analysis_3d.py
index 23cbda24a13..4f9e75e11d7 100755
--- a/Examples/Tests/langmuir/analysis_3d.py
+++ b/Examples/Tests/langmuir/analysis_3d.py
@@ -25,6 +25,7 @@
 
 import numpy as np
 from analysis_utils import check_charge_conservation
+from openpmd_viewer import OpenPMDTimeSeries
 from scipy.constants import c, e, epsilon_0, m_e
 
 # test name
@@ -57,15 +58,19 @@
 cos = {"Ex": (0, 1, 1), "Ey": (1, 0, 1), "Ez": (1, 1, 0)}
 
 
-def get_contribution(is_cos, k, idim):
-    du = (hi[idim] - lo[idim]) / Ncell[idim]
-    u = lo[idim] + du * (0.5 + np.arange(Ncell[idim]))
+def get_contribution_at_positions(is_cos, k, idim, u):
     if is_cos[idim] == 1:
         return np.cos(k * u)
     else:
         return np.sin(k * u)
 
 
+def get_contribution(is_cos, k, idim):
+    du = (hi[idim] - lo[idim]) / Ncell[idim]
+    u = lo[idim] + du * (0.5 + np.arange(Ncell[idim]))
+    return get_contribution_at_positions(is_cos, k, idim, u)
+
+
 def get_theoretical_field(field, t):
     amplitude = epsilon * (m_e * c**2 * k[field]) / e * np.sin(wp * t)
     cos_flag = cos[field]
@@ -83,6 +88,18 @@ def get_theoretical_field(field, t):
     return E
 
 
+def get_theoretical_field_at_positions(field, t, x, y, z):
+    amplitude = epsilon * (m_e * c**2 * k[field]) / e * np.sin(wp * t)
+    cos_flag = cos[field]
+    x_contribution = get_contribution_at_positions(cos_flag, kx, 0, x)
+    y_contribution = get_contribution_at_positions(cos_flag, ky, 1, y)
+    z_contribution = get_contribution_at_positions(cos_flag, kz, 2, z)
+
+    E = amplitude * x_contribution * y_contribution * z_contribution
+
+    return E
+
+
 # Read the file
 ds = yt.load(fn)
 
@@ -122,6 +139,17 @@ def get_theoretical_field(field, t):
     print("%s: Max error: %.2e" % (field, max_error))
     error_rel = max(error_rel, max_error)
 
+ts = OpenPMDTimeSeries("./diags/openpmd")
+x, y, z = ts.get_particle(["x", "y", "z"], species="electrons", iteration=40)
+for field in ["Ex", "Ey", "Ez"]:
+    E_sim_particles = ts.get_particle(
+        [field.lower()], species="electrons", iteration=40
+    )
+    E_th_particles = get_theoretical_field_at_positions(field, t0, x, y, z)
+    max_error = abs(E_sim_particles - E_th_particles).max() / abs(E_th_particles).max()
+    print("%s: Max error at particles: %.2e" % (field, max_error))
+    error_rel = max(error_rel, max_error)
+
 # Plot the last field from the loop (Ez at iteration 40)
 fig, (ax1, ax2) = plt.subplots(1, 2, dpi=100)
 # First plot (slice at y=0)
diff --git a/Examples/Tests/langmuir/inputs_base_3d b/Examples/Tests/langmuir/inputs_base_3d
index 9b8b9b16e41..40ec3852f65 100644
--- a/Examples/Tests/langmuir/inputs_base_3d
+++ b/Examples/Tests/langmuir/inputs_base_3d
@@ -90,9 +90,14 @@ positrons.momentum_function_uz(x,y,z) = "-epsilon * k/kp * cos(k*x) * cos(k*y) *
 
 # Diagnostics
 warpx.synchronize_velocity_for_diagnostics = 1
-diagnostics.diags_names = diag1
+diagnostics.diags_names = diag1 openpmd
 diag1.intervals = max_step
 diag1.diag_type = Full
 diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell rho divE
 diag1.electrons.variables = x y z w ux
 diag1.positrons.variables = x y z uz
+
+openpmd.intervals = 40
+openpmd.diag_type = Full
+openpmd.format = openpmd
+openpmd.electrons.additional_variables = Ex Ey Ez
diff --git a/Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_JRhom_LL2_picmi.py b/Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_JRhom_LL2_picmi.py
index 76caccf7d62..acdc47d04d8 100755
--- a/Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_JRhom_LL2_picmi.py
+++ b/Examples/Tests/langmuir/inputs_test_3d_langmuir_multi_psatd_JRhom_LL2_picmi.py
@@ -135,8 +135,15 @@
     species=[positrons],
     data_list=["x", "y", "z", "uz"],
 )
+field_on_particle_diag = picmi.ParticleDiagnostic(
+    name="openpmd",
+    period=40,
+    data_list=["x", "y", "z", "Ex", "Ey", "Ez"],
+    warpx_format="openpmd",
+)
 sim.add_diagnostic(particle_diag1_electrons)
 sim.add_diagnostic(particle_diag1_positrons)
+sim.add_diagnostic(field_on_particle_diag)
 
 sim.write_input_file(file_name="inputs_test_3d_langmuir_multi_psatd_JRhom_LL2_picmi")
 
diff --git a/Regression/Checksum/benchmarks_json/test_2d_particle_EM_diagnostics.json b/Regression/Checksum/benchmarks_json/test_2d_particle_EM_diagnostics.json
new file mode 100644
index 00000000000..d4f1470f835
--- /dev/null
+++ b/Regression/Checksum/benchmarks_json/test_2d_particle_EM_diagnostics.json
@@ -0,0 +1,25 @@
+{
+  "lev=0": {
+    "Bx": 3864.6903326868182,
+    "By": 239757.36697092294,
+    "Bz": 119878.42982987872,
+    "Ex": 2396689176285.2944,
+    "Ey": 35934018061853.67,
+    "Ez": 71884648325011.86
+  },
+  "electron": {
+    "particle_bx": 3.890980338189904e-07,
+    "particle_by": 0.34760630481904997,
+    "particle_bz": 0.17371926374000385,
+    "particle_ex": 0.0001349025551462546,
+    "particle_ey": 50608854.55959178,
+    "particle_ez": 111088826.05245684,
+    "particle_position_x": 2.5e-06,
+    "particle_position_y": 0.0,
+    "particle_position_z": 1.25e-05,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_weight": 0.0
+  }
+}
diff --git a/Source/Diagnostics/ParticleDiag/ParticleDiag.H b/Source/Diagnostics/ParticleDiag/ParticleDiag.H
index 18b6e450048..7a0bd871bef 100644
--- a/Source/Diagnostics/ParticleDiag/ParticleDiag.H
+++ b/Source/Diagnostics/ParticleDiag/ParticleDiag.H
@@ -24,6 +24,12 @@ public:
     [[nodiscard]] std::string getSpeciesName() const { return m_name; }
     amrex::Vector<int> m_plot_flags;
     bool m_plot_phi = false; // Whether to output the potential phi on the particles
+    bool m_plot_Ex = false; // Whether to output the Ex on the particles
+    bool m_plot_Ey = false; // Whether to output the Ey on the particles
+    bool m_plot_Ez = false; // Whether to output the Ez on the particles
+    bool m_plot_Bx = false; // Whether to output the Bx on the particles
+    bool m_plot_By = false; // Whether to output the By on the particles
+    bool m_plot_Bz = false; // Whether to output the Bz on the particles
 
     bool m_do_random_filter  = false;
     bool m_do_uniform_filter = false;
diff --git a/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp b/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp
index e90278d9cbe..8a72e1015d7 100644
--- a/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp
+++ b/Source/Diagnostics/ParticleDiag/ParticleDiag.cpp
@@ -46,11 +46,15 @@ ParticleDiag::ParticleDiag (
                 if (var == "z") { var = "phi"; }
 #endif
                 if (var == "phi") {
-                    // User requests phi on particle. This is *not* part of the variables that
-                    // the particle container carries, and is only added to particles during output.
-                    // Therefore, this case needs to be treated specifically.
                     m_plot_phi = true;
-                } else {
+                }
+                else if (var == "Ex") { m_plot_Ex = true; }
+                else if (var == "Ey") { m_plot_Ey = true; }
+                else if (var == "Ez") { m_plot_Ez = true; }
+                else if (var == "Bx") { m_plot_Bx = true; }
+                else if (var == "By") { m_plot_By = true; }
+                else if (var == "Bz") { m_plot_Bz = true; }
+                else {
                     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(pc->HasRealComp(var),
                         "variables argument '" + var
                         +"' is not an existing attribute for this species");
@@ -72,6 +76,22 @@ ParticleDiag::ParticleDiag (
         }
     }
 
+    amrex::Vector<std::string> additional_variables;
+    pp_diag_name_species_name.queryarr("additional_variables", additional_variables);
+    for (auto& var : additional_variables) {
+        if (var == "phi") { m_plot_phi = true; }
+        else if (var == "Ex") { m_plot_Ex = true; }
+        else if (var == "Ey") { m_plot_Ey = true; }
+        else if (var == "Ez") { m_plot_Ez = true; }
+        else if (var == "Bx") { m_plot_Bx = true; }
+        else if (var == "By") { m_plot_By = true; }
+        else if (var == "Bz") { m_plot_Bz = true; }
+        else {
+            WARPX_ABORT_WITH_MESSAGE(
+                "additional_variables argument '" + var + "' is not a valid variable");
+        }
+    }
+
 #if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
     // Always write out theta, whether or not it's requested,
     // to be consistent with always writing out r and z.
diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp
index ac19f98459c..e48e67b7042 100644
--- a/Source/Diagnostics/ParticleIO.cpp
+++ b/Source/Diagnostics/ParticleIO.cpp
@@ -7,12 +7,13 @@
  * License: BSD-3-Clause-LBNL
  */
 
-#include "Fields.H"
 #include "Particles/ParticleIO.H"
+
+#include "Fields.H"
+#include "Particles/LaserParticleContainer.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 #include "Particles/MultiParticleContainer.H"
 #include "Particles/PhysicalParticleContainer.H"
-#include "Particles/LaserParticleContainer.H"
 #include "Particles/RigidInjectedParticleContainer.H"
 #include "Particles/WarpXParticleContainer.H"
 #include "Utils/TextMsg.H"
@@ -303,3 +304,121 @@ storePhiOnParticles ( WarpXParticleContainer::Base& tmp,
         }
     }
 }
+
+void
+storeFieldOnParticles ( WarpXParticleContainer::Base& tmp, bool is_full_diagnostic,
+                        bool save_Ex, bool save_Ey, bool save_Ez,
+                        bool save_Bx, bool save_By, bool save_Bz) {
+
+    using PinnedParIter = typename WarpXParticleContainer::ParIterType;
+    using ablastr::fields::Direction;
+
+    // When this is not a full diagnostic, the particles are not written at the same physical time (i.e. PIC iteration)
+    // that they were collected. This happens for diagnostics that use buffering (e.g. BackTransformed, BoundaryScraping).
+    // Here the field is gathered at the iteration when particles are written (not collected) and is thus mismatched.
+    // To avoid confusion, we raise an error in this case.
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        is_full_diagnostic,
+        "Output of the fields on the particles was requested, "
+        "but this is only available with `diag_type = Full`.");
+
+    if (save_Ex) { tmp.AddRealComp("Ex"); }
+    if (save_Ey) { tmp.AddRealComp("Ey"); }
+    if (save_Ez) { tmp.AddRealComp("Ez"); }
+    if (save_Bx) { tmp.AddRealComp("Bx"); }
+    if (save_By) { tmp.AddRealComp("By"); }
+    if (save_Bz) { tmp.AddRealComp("Bz"); }
+
+    int const Ex_index = save_Ex ? tmp.GetRealCompIndex("Ex") : -1;
+    int const Ey_index = save_Ey ? tmp.GetRealCompIndex("Ey") : -1;
+    int const Ez_index = save_Ez ? tmp.GetRealCompIndex("Ez") : -1;
+    int const Bx_index = save_Bx ? tmp.GetRealCompIndex("Bx") : -1;
+    int const By_index = save_By ? tmp.GetRealCompIndex("By") : -1;
+    int const Bz_index = save_Bz ? tmp.GetRealCompIndex("Bz") : -1;
+
+    auto & warpx = WarpX::GetInstance();
+    auto & fields = warpx.m_fields;
+
+    const int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes;
+    const int nox = WarpX::nox;
+    const bool galerkin_interpolation = WarpX::galerkin_interpolation;
+
+    for (int lev=0; lev<=warpx.finestLevel(); lev++) {
+
+        const amrex::XDim3 dinv = WarpX::InvCellSize(lev);
+
+        amrex::MultiFab & Ex = *fields.get(FieldType::Efield_aux, Direction{0}, lev);
+        amrex::MultiFab & Ey = *fields.get(FieldType::Efield_aux, Direction{1}, lev);
+        amrex::MultiFab & Ez = *fields.get(FieldType::Efield_aux, Direction{2}, lev);
+        amrex::MultiFab & Bx = *fields.get(FieldType::Bfield_aux, Direction{0}, lev);
+        amrex::MultiFab & By = *fields.get(FieldType::Bfield_aux, Direction{1}, lev);
+        amrex::MultiFab & Bz = *fields.get(FieldType::Bfield_aux, Direction{2}, lev);
+
+#ifdef AMREX_USE_OMP
+        #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
+#endif
+        for (PinnedParIter pti(tmp, lev); pti.isValid(); ++pti) {
+
+            amrex::Box box = pti.tilebox();
+            const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, lev, 0._rt);
+            const Dim3 lo = lbound(box);
+
+            FArrayBox const* exfab = &Ex[pti];
+            FArrayBox const* eyfab = &Ey[pti];
+            FArrayBox const* ezfab = &Ez[pti];
+            FArrayBox const* bxfab = &Bx[pti];
+            FArrayBox const* byfab = &By[pti];
+            FArrayBox const* bzfab = &Bz[pti];
+
+            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();
+
+            const auto getPosition = GetParticlePosition<PIdx>(pti);
+
+            amrex::ParticleReal* Ex_particle_arr = save_Ex ? pti.GetStructOfArrays().GetRealData(Ex_index).dataPtr() : nullptr;
+            amrex::ParticleReal* Ey_particle_arr = save_Ey ? pti.GetStructOfArrays().GetRealData(Ey_index).dataPtr() : nullptr;
+            amrex::ParticleReal* Ez_particle_arr = save_Ez ? pti.GetStructOfArrays().GetRealData(Ez_index).dataPtr() : nullptr;
+            amrex::ParticleReal* Bx_particle_arr = save_Bx ? pti.GetStructOfArrays().GetRealData(Bx_index).dataPtr() : nullptr;
+            amrex::ParticleReal* By_particle_arr = save_By ? pti.GetStructOfArrays().GetRealData(By_index).dataPtr() : nullptr;
+            amrex::ParticleReal* Bz_particle_arr = save_Bz ? pti.GetStructOfArrays().GetRealData(Bz_index).dataPtr() : nullptr;
+
+            // Loop over the particles and update their position
+            amrex::ParallelFor( pti.numParticles(),
+                [=] AMREX_GPU_DEVICE (long ip) {
+
+                    amrex::ParticleReal xp, yp, zp;
+                    getPosition(ip, xp, yp, zp);
+
+                    amrex::ParticleReal Exp = 0., Eyp = 0., Ezp = 0.;
+                    amrex::ParticleReal Bxp = 0., Byp = 0., Bzp = 0.;
+
+                    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,
+                                   dinv, xyzmin, lo, n_rz_azimuthal_modes,
+                                   nox, galerkin_interpolation);
+
+
+                    if (Ex_particle_arr) Ex_particle_arr[ip] = Exp;
+                    if (Ey_particle_arr) Ey_particle_arr[ip] = Eyp;
+                    if (Ez_particle_arr) Ez_particle_arr[ip] = Ezp;
+                    if (Bx_particle_arr) Bx_particle_arr[ip] = Bxp;
+                    if (By_particle_arr) By_particle_arr[ip] = Byp;
+                    if (Bz_particle_arr) Bz_particle_arr[ip] = Bzp;
+
+                }
+            );
+        }
+    }
+}
diff --git a/Source/Diagnostics/WarpXOpenPMD.cpp b/Source/Diagnostics/WarpXOpenPMD.cpp
index 97522101ac6..52070cbdeb4 100644
--- a/Source/Diagnostics/WarpXOpenPMD.cpp
+++ b/Source/Diagnostics/WarpXOpenPMD.cpp
@@ -605,8 +605,17 @@ for (const auto & particle_diag : particle_diags) {
     }
 
     // Gather the electrostatic potential (phi) on the macroparticles
+    // Detect whether this is a FullDiagnostic (use_pinned_pc = 0) or
+    // a BoundaryScrapingDiagnostic/BackTransformedDiagnostic (use_pinned_pc=1)
+    bool const is_full_diagnostic = !use_pinned_pc;
     if ( particle_diag.m_plot_phi ) {
-        storePhiOnParticles( tmp, WarpX::electrostatic_solver_id, !use_pinned_pc );
+        storePhiOnParticles( tmp, WarpX::electrostatic_solver_id, is_full_diagnostic);
+    }
+    if (particle_diag.m_plot_Ex || particle_diag.m_plot_Ey || particle_diag.m_plot_Ez ||
+        particle_diag.m_plot_Bx || particle_diag.m_plot_By || particle_diag.m_plot_Bz) {
+        storeFieldOnParticles(tmp, is_full_diagnostic,
+                              particle_diag.m_plot_Ex, particle_diag.m_plot_Ey, particle_diag.m_plot_Ez,
+                              particle_diag.m_plot_Bx, particle_diag.m_plot_By, particle_diag.m_plot_Bz);
     }
 
     // names of amrex::ParticleReal and int particle attributes in SoA data
diff --git a/Source/Particles/ParticleIO.H b/Source/Particles/ParticleIO.H
index 72b7d363647..3bfd8a88d1b 100644
--- a/Source/Particles/ParticleIO.H
+++ b/Source/Particles/ParticleIO.H
@@ -91,4 +91,21 @@ void
 storePhiOnParticles ( WarpXParticleContainer::Base& tmp,
     ElectrostaticSolverAlgo electrostatic_solver_id, bool is_full_diagnostic );
 
+/** \brief Gathers fields from a MultiFab to the macroparticles.
+ *         Adds a runtime component of the particle container to store it.
+ *
+ * \param tmp The particle container on which to store the gathered field
+ * \param is_full_diagnostic Whether this diagnostic is a full diagnostic
+ * \param save_Ex Whether the Ex field is saved
+ * \param save_Ey Whether the Ey field is saved
+ * \param save_Ez Whether the Ez field is saved
+ * \param save_Bx Whether the Bx field is saved
+ * \param save_By Whether the By field is saved
+ * \param save_Bz Whether the Bz field is saved
+ */
+void
+storeFieldOnParticles (WarpXParticleContainer::Base& tmp, bool is_full_diagnostic,
+                       bool save_Ex, bool save_Ey, bool save_Ez,
+                       bool save_Bx, bool save_By, bool save_Bz);
+
 #endif /* WARPX_PARTICLEIO_H_ */