Another thing you might try is Nearest.

Nearest[deltaSequence -> {"Index", "Distance"}, particleValues, {All, deltaParam}]

I renamed some variables because I was having trouble keeping track, but I think it would be this with your original variables:

Nearest[xg -> {"Index", "Distance"}, xp, {All, dx}]

I included my renamed version so that if I messed it up, you can hopefully correct it.

Anyway, what this will give you is something like the following:

{{{6, 0.0122529}, {7, 0.0877471}}, {{18, 0.00370914}, {17, 
   0.0962909}}, {{18, 0.0247811}, {17, 0.0752189}}, {{4, 
   0.0364263}, {5, 0.0635737}}, {{8, 0.0122352}, {7, 0.0877648}}, <...etc...>}

You should be able to write a function that uses this data to do your update to your rho table. I haven't done a performance comparison.

Related question on MSE.

I am confused. so you are saying my code should be something like this?

Do[pos=Floor[x[[i]]/dx];rho[[pos]]=q/dx * Abs[FractionalPart[x[[i]]/dx], {i,nparticles}];
 Do[pos=Ceiling[x[[i]]/dx];\[rho][[pos]]=q/dx * Abs[FractionalPart[x[[i]]/dx], {i,nparticles}];

Because the particles influence the point of the grid before and after them.

Yes. ρ is calculated at each point xg[[j]]. The number of times ρ is updated depends on the number of particles that has a position between xg[[j]] + or - dx.