I’m working basically with the creation of random numbers that represent the positions of $N$ electrons within a box. Later, I sum all the individual electric fields generated by the electrons, thus yielding the total electric field at one specific point in space (outside the box). After doing that for $N_{re}$ realizations of random electrons position, I build a histogram for the electric field at that specific point.

When I chose $N=10000$ electrons and $N_{re}=10000$, my computer takes a lot of time to perform these calculations — even using Parallel commands.

Would you guys have any suggestion for improvement on the speed calculation? I’d appreciate anything that might help.

Here is the code, it’s actually very short. I should notice that I’m also fitting the histogram with a well known distribution.

```
MonteCarlo(k_, Nt_, R_, Nre_, z_) := Module({el, (Epsilon)d, nI, Rdip},
(Epsilon)d = 9.6 8.8542 10^-12;
el = 1.6 10^-19;
nI = Nt/R^3;
rdop(i_) := {RandomReal({-R/2, R/2}, {1})((1)),RandomReal({-R/2, R/2}, {1})((1)),RandomReal({-R/2, R/2}, {1})((1))};
Rlistdop = ParallelTable(rdop(i), {i, 1, Nt});
AbsoluteTiming(Do(Rlistdip(j) = ParallelTable(rdip(i), {i, 1, Nt}), {j, 1,Nre}));
El(j_, i_):=el/(4 (Pi) (Epsilon)d) (Rlistdop((i)) - rNV)/(Norm(Rlistdop((i))-rNV)^3);
rNV = {0, 0, -R/2 - z};
Et(j_) := Sum(El(j, i), {i, 1, Nt});
Elmean(j_) := 1/Nt Sum(El(j, i), {i, 1, Nt});
dataEz(k)=ParallelTable(If(Abs(Et(j)((3))) < 1 10^8, Et(j)((3)), Nothing), {j, 1, Nre});
pEz(k)=Histogram({dataEz(k)}, {-2 10^7, 2 10^7, 10^4}, PDF,PlotRange -> {{-2 10^7, 2 10^7},
All}, PlotTheme -> "Scientific",ChartLayout -> "Overlapped",LabelStyle ->Directive(Black,
{FontFamily -> "Latin Modern Roman", FontSize -> 17}));
(Mu)0(k)=(Mu)(k)/.FindDistributionParameters(dataEz(k),StudentTDistribution((Mu)(k),bba(k),cca(k)))((1));
b0(k)=bba(k)/.FindDistributionParameters(dataEz(k),StudentTDistribution((Mu)(k),bba(k),cca(k)))((2));
c0(k)=cca(k)/.FindDistributionParameters(dataEz(k),StudentTDistribution((Mu)(k),bba(k),cca(k)))((3));
(ScriptCapitalD)(k) = StudentTDistribution((Mu)0(k), b0(k), c0(k));
(Mu)0(k) = Median((ScriptCapitalD)(k));
fitd(x_)(k) := PDF((ScriptCapitalD)(k), x);
xhalf1(k)=x/.Solve(PDF((ScriptCapitalD)(k), (Mu)0(k))/2 ==PDF((ScriptCapitalD)(k),x),x)((1));
xhalf2(k)=x/.Solve(PDF((ScriptCapitalD)(k),(Mu)0(k))/2==PDF((ScriptCapitalD)(k),x),x)((2));
Ehalf(k) = Abs(xhalf1(k)) + Abs(xhalf2(k));
Ehalfa(k, z) = Ehalf(k);
)
```

Using the module above, my situation correspond on

```
Nt = 10000;
R = 1. 10^-7;
Nre = 10000;
MonteCarlo(1, Nt, R, Nre, 2 R)
```

Thanks!

Best,

Denis