Leapfrog method in Particle-in-cell Recently, I wanted to write a 3D electromagnetic Particle-in-cell code with C++. I know that I should use leapfrog method, for example, when I calculate the position and velocity of particles，I should calculate according to the following order:

//firt move velocity back by 0.5dt
v = v - 0.5*q/m*E*dt;
//main loop
//update velocity and position
v = v + q/m*E*dt;
x = x + v*dt;

But when I want to write 3D PIC, I don't know the calculation order of $E,\,B,\,v,\,x$ (postion) and their formula. Can someone help me please?
 A: For any PIC simulation, you are necessarily tying yourself to particles, particles that experience forces. Thus, we have the generic force law:
\begin{align}
m_i\frac{d\mathbf v_i}{dt}&=\mathbf F_i  \tag{1a}\\
\frac{d\mathbf x_i}{dt}&=\mathbf v_i \tag{1b}
\end{align}
In the case of PIC, you are often considering the electromagnetic (Lorentz) force, so the force on particle $i$ is,
$$
\mathbf F_i=q_i\left(\mathbf E+\boldsymbol\beta_i\times\mathbf B\right)\tag{2}
$$
with $\boldsymbol\beta=\mathbf v/c$. Since we've got particles, the electric and magnetic fields are defined as
$$
\mathbf E=\sum_{i\neq j}\frac{q_j}{r_{ij}^2}\hat{r}_j,\qquad\mathbf B=\sum_{j\neq i}\boldsymbol\beta_j\times\mathbf E
$$
You've already mentioned the leapfrog method, which is one of the more common methods for PIC simulations (at least to my knowledge). Your question is how/when to compute $\mathbf E,\,\mathbf B$. The answer is quite simple: after you've computed $x$; after all, the electric and magnetic fields are functions of space.
Note, though, that the magnetic field term of the Lorentz force (eq (2)) contains a velocity component; this must also be leapfrogged:
$$
\mathbf v^{n+1/2}_i\sim\frac12\left(\mathbf v_i^{n-1/2}+\mathbf v_i^{n+1/2}\right)\times\mathbf B
$$
which that $\mathbf v_i^{n+1/2}$ on the RHS makes it a little bit harder to solve (you'll likely need a linear algebra solver). There is an alternative, called the Boris method, that combines the velocity and electric field into a single term & you end up with a "rotation" without needing to involve the linear algebra solver needed; however, I am not well versed in this so I cannot say much beyond that.
Thus, the generic stepping algorithm would be something like,
while (t < tend)
   computeEMFields(x, v, E, B);
   newV = updateVelocity(v, E, B);
   newX = moveParticles(v, newV, x);
   computeDT();
   t = t + dt;
   x = newX
   v = newV
end while

