# Monte Carlo Simulation of Charged Particles in Non-Uniform Electric Field

I have a code (provided below) which simulates the motion of an ensemble of charged particles which are subjected to a complex static electric field which passes through some aperture in a metallic surface (N spheres, whose number is defined as numbodies in my code).

Now, I would like to run the simulation for a single particle many, many times, and plot each of their final resting places on the surface (z=-partRad) of my metalized aperture at (x(tfin), y(tfin), z= -partRad) together, in 2D. These coordinates can be displayed as simple dots in the x-y plane.

It would also be cool to superimpose the surface of my aperture onto the x-y graph of the final particle locations, to see exactly where they most frequently land. I would also like to find the average distance these particles center from some user specified point(s), after everything is said and done. I would also appreciate any cool ideas for statistical analysis regarding the outcome of the Monte Carlo simulation!!!

Each of the particles should be given a random starting position and velocity given by

vel0 = Table[
Partition[{RandomReal[{-0.0001, 0.0001}],
RandomReal[{-0.0001, 0.0001}], RandomReal[{-0.0001, 0.0001}]}, 1],
numbodies]
pos0 = Table[
Partition[{RandomReal[{0.5, 1}], RandomReal[{0.5, 1}],
RandomReal[{-0.4, -0.04}]}, 1], numbodies];


in my code.

The first part of the code, above the (***... MONTE CARLO CAN START...***) line in my code, doesn't need to be included in the Monte Carlo simulation (unless you find some reason I do not see).

Here is the code. Happy hunting!

Clear["Global*"];
Needs["NDSolveFEM"]

q = 1.617733*10^-18;(*Net ion Charge*)

R = Import[
"https://www.dropbox.com/s/dds8rm3odg2m7gu/largeAp.obj?dl=1"];
RegionDimension[R];
M = BoundaryMeshRegion[MeshCoordinates[R], MeshCells[R, 2]];
RegionDimension[M];
Volume[M];

r = RegionDifference[
RegionDifference[
RegionDifference[Cuboid[{0, 0, -6.5}, {2, 2, 0.5}], M],
Cuboid[{0, 0, 0.4}, {2, 2, 0.5}]],
Cuboid[{0, 0, -6.5}, {2, 2, -6.4}]];
ToElementMesh[r]["Wireframe"];
pol = -1;

V0 = 1000;

sol = NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] == -pol*
V0/2, (0.4 <= z <= 0.5) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
pol*V0/2, (0.0071 <= z <= 0.0072) && (0 <= y <= 2) && (0 <= x <=
2)], DirichletCondition[
V[x, y, z] ==
0, (0 <= z <= 0.0070) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
0, (-6.5 <= z <= -6.4) && (0 <= y <= 2) && (0 <= x <= 2)]},
V, {x, y, z} \[Element] r];

electricField[x_, y_, z_] = -Grad[sol[x, y, z], {x, y, z}];

v = Show[VectorPlot3D[
electricField[x, y, z], {x, 0.5, 1}, {y, 0.5, 1}, {z, -0.5, 0.1},
PlotTheme -> "Detailed", ColorFunction -> "Rainbow",
PerformanceGoal -> "Quality", VectorScale -> 0.05,
VectorPoints -> 7, PlotLegends -> Automatic], M];

mass = 6.52*10^-11;(*particle mass in kg/m^3*)

eforceX[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][];
eforceY[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][];
eforceZ[x_, y_, z_] = q*Grad[sol[x, y, z], {x, y, z}][];

(***********************    MONTE CARLO CAN START    ******************************)

numbodies = 1;

vel0 = Table[
Partition[{RandomReal[{-0.0001, 0.0001}],
RandomReal[{-0.0001, 0.0001}], RandomReal[{-0.0001, 0.0001}]}, 1],
numbodies]
pos0 = Table[
Partition[{RandomReal[{0.5, 1}], RandomReal[{0.5, 1}],
RandomReal[{-0.4, -0.04}]}, 1], numbodies];

tfin = 1000;

eqs = Table[{x1[j]''[t] == -a[t]*(1/mass)*eforceX[x, y, z] /. {x ->
x1[j][t], y -> y1[j][t], z -> z1[j][t]},
y1[j]''[t] == -a[t]*(1/mass)*eforceY[x, y, z] /. {x -> x1[j][t],
y -> y1[j][t], z -> z1[j][t]},
z1[j]''[t] == -a[t]*(1/mass)*eforceZ[x, y, z] /. {x -> x1[j][t],
y -> y1[j][t], z -> z1[j][t]}, x1[j] == pos0[[j, 1, 1]],
y1[j] == pos0[[j, 2, 1]], z1[j] == pos0[[j, 3, 1]],
x1[j]' == vel0[[j, 1, 1]], y1[j]' == vel0[[j, 2, 1]],
z1[j]' == vel0[[j, 3, 1]], a == 1}, {j, numbodies}];

vars = Flatten[Table[{x1[j][t], y1[j][t], z1[j][t]}, {j, numbodies}]];

event = Table[
WhenEvent[
z1[j][t] >= -partRad, {a[t] -> 0, x1[j]'[t] -> 0, y1[j]'[t] -> 0,
z1[j]'[t] -> 0}] /. {j -> i}, {i, numbodies}];

sol1 = NDSolve[{eqs, event}, Head /@ vars, {t, 0, tfin},
DiscreteVariables -> a]

frames = Table[
Show[v, ParametricPlot3D[
Table[{x1[j][t], y1[j][t], z1[j][t]} /. sol1, {j,
numbodies}], {t, 0, tf}, PlotRange -> All, Axes -> Off],
Graphics3D[
Table[{Hue[.35],
Sphere[{x1[j][tf], y1[j][tf], z1[j][tf]} /. sol1, 0.025]}, {j,
numbodies}]]], {tf, .015*tfin, tfin, .015*tfin}];

ListAnimate[Rasterize[#, "Image"] & /@ frames]

• Hi @DanielHuber. Are you running into a problem trying to execute the code? I copy/pasted the exact code provided above, and it runs without error on my end.
– Zach
Aug 1 at 14:02
• Sorry, I must have done something wrong by copy/past. There was a Table command without iterator. Aug 1 at 15:16
• @zak720 As I remember we have discussed once this mesh with two holes on mathematica.stackexchange.com/questions/269806/… Aug 2 at 10:03

To simulate hundreds and thousands runs we need very fast code, while at present form code for trajectories is very slow. To improve code, first, we include the electric field computation in NDSolve as follows

Clear["Global*"];
Needs["NDSolveFEM"]

q = 1.617733*10^-18;(*Net ion Charge*)

R = Import[
"https://www.dropbox.com/s/dds8rm3odg2m7gu/largeAp.obj?dl=1"];
RegionDimension[R];
M = BoundaryMeshRegion[MeshCoordinates[R], MeshCells[R, 2]];
r = RegionDifference[
RegionDifference[
RegionDifference[Cuboid[{0, 0, -6.5}, {2, 2, 0.5}], M],
Cuboid[{0, 0, 0.4}, {2, 2, 0.5}]],
Cuboid[{0, 0, -6.5}, {2, 2, -6.4}]];
pol = -1;

V0 = 1000;

{U, {Fx, Fy, Fz}} =
NDSolveValue[{Laplacian[V[x, y, z], {x, y, z}] == 0,
DirichletCondition[
V[x, y, z] == -pol*
V0/2, (0.4 <= z <= 0.5) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
pol*V0/2, (0.0071 <= z <= 0.0072) && (0 <= y <= 2) && (0 <= x <=
2)], DirichletCondition[
V[x, y, z] ==
0, (0 <= z <= 0.0070) && (0 <= y <= 2) && (0 <= x <= 2)],
DirichletCondition[
V[x, y, z] ==
0, (-6.5 <= z <= -6.4) && (0 <= y <= 2) && (0 <= x <=
2)]}, {V, -Grad[V[x, y, z], {x, y, z}]}, {x, y, z} \[Element]
r];


Visualization

v = Show[
VectorPlot3D[{Fx, Fy, Fz}, {x, 0.5, 1}, {y, 0.5, 1}, {z, -0.5, 0.1},
PlotTheme -> "Detailed", ColorFunction -> "Rainbow",
PerformanceGoal -> "Quality", VectorScale -> 0.05,
VectorPoints -> 7, PlotLegends -> Automatic],
RegionPlot3D[R, PlotStyle -> Opacity[0.5]]] Second, we use Verlet integration algorithm to simulate trajectories, we have

mass = 6.52*10^-11;(*particle mass*)
SeedRandom;
n = 100; X0 = Take[RandomReal[{.5, 1}, 2 n], n]; Y0 =
Take[RandomReal[{.5, 1}, 2 n], -n]; Z0 =
RandomReal[{-.4, -.04}, n]; rrv =
RandomReal[{-0.0001, 0.0001}, 3 n]; u0 = Take[rrv, n]; v0 =
Take[rrv, {n + 1, 2 n}]; w0 = Take[rrv, -n];

SymplecticLeapfrog = {"SymplecticPartitionedRungeKutta",
"DifferenceOrder" -> 2, "PositionVariables" :> qvars}; time = {t, 0,
2}; tm = 10^3; partRad = 0.0125; qvars = {X[t], Y[t], Z[t]};
eqn[i_] := {X''[t] == q/mass tm^2 Fx /. {x -> X[t], y -> Y[t],
z -> Z[t]},
Y''[t] == q/mass tm^2 Fy /. {x -> X[t], y -> Y[t], z -> Z[t]},
Z''[t] == q/mass tm^2 Fz /. {x -> X[t], y -> Y[t], z -> Z[t]},
X == X0[[i]], Y == Y0[[i]], Z == Z0[[i]], X' == u0[[i]],
Y' == v0[[i]], Z' == w0[[i]],
WhenEvent[Z[t] >= -partRad, end = t; "StopIntegration"]}; h = 10^-4;
Do[soln[i] =
NDSolve[eqn[i], qvars, time, StartingStepSize -> h,
Method -> SymplecticLeapfrog];
T[i] = end;, {i, 1, n}]


It takes about 25 sec to compute 100 trajectories on my laptop. Visualization of trajectories

Show[v, Table[
ParametricPlot3D[qvars /. soln[i], {t, 0, T[i]},
ColorFunction -> Hue], {i, 100}]] Finally we can compute point distribution as

fin = Table[qvars /. soln[i][] /. t -> T[i], {i, n}];
Show[Graphics3D[{{Red, Point[fin]}},
PlotRange -> {{.5, 1}, {.5, 1}, All}],
RegionPlot3D[R, PlotStyle -> Opacity[0.5]]] To compute statistics we need 10^4 points as minimum. Example with 1000 points Example with $$10^4$$ points In the last case 33 particles of 10000 are not deposited on the plane z==-partRad, therefore we select cases and plot as follows

fin = Table[qvars /. soln[i][] /. t -> T[i], {i, n}];

fin1 = Select[fin, Norm[#[] + partRad] <= .01 &];

Show[Graphics3D[{{Red, Point[fin1]}},
PlotRange -> {{.5, 1}, {.5, 1}, All}],
RegionPlot3D[R, PlotStyle -> Opacity[0.5]]]


Using fin1 we can compute mean and standard deviation on every coordinate

XYmean = Mean[fin1]

(*Out[]= {0.757097, 0.754464, -0.0125015}*)

std = StandardDeviation[fin1]

(*Out[]= {0.0511615, 0.045477, 0.000108942}*)


Visualization of mean and std[] as a sphere

Show[Graphics3D[{{Red, Point[fin1]}, {Blue, Opacity[0.25],
Sphere[XYmean, std[]]}}, PlotRange -> {{.5, 1}, {.5, 1}, All},
Axes -> {True, True, False}],
RegionPlot3D[R, PlotStyle -> Opacity[0.5]]] • this is an absolute work of art. I had come up with a very crude version of this prior to seeing this answer (and it is, as you might guess, extremely slow). Fantastic work.
– Zach
Aug 5 at 1:11
• How would I set a WhenEvent such that when a particle comes into contact with any other particle position previously simulated, its motion comes to a stop? In my crude version, I made an empty list called outCoord and filled it with the final resting coordinates of each particle simulated. I then check whether, WhenEvent[ Sqrt[(outCoord[[k, 1]][[1, 1]] - x1[j][t])^2 + (outCoord[[k, 1]][[1, 2]] - y1[j][t])^2 + (outCoord[[k, 1]][[1, 3]] - z1[j][t])^2] <= partRad for each new particle simulated, such that if this is true, it stops. How can I do it here?
– Zach
Aug 5 at 3:19
• I don't understand last question. Do you mean collisions or something else? Aug 5 at 7:56
• Yes like a collision. If a particle encounters a region of space where previously simulated particles have landed, it's motion stops. This would create a 3D structure of particles on the surface since they would now be able to "stack" on top of eachother.
– Zach
Aug 5 at 12:37
• Ok, for this we need to define some size of particle and rule for packaging. Aug 5 at 12:43