# Calculating minimum potential energies for Thomson problem

I have found a table of values for minimum potential energies up to 470 electrons on Wikipedia. However, the method I used for calculating energies is not efficient enough for configurations that have more than 20 electrons.

Could someone optimize the code or provide some insights for optimization, please? (I only need precision up to 2 decimal places)

Clear["Global*"]
force[pt1_, pt2_] :=
If[pt1 == pt2, {0, 0, 0},
Normalize[(pt1 - pt2)]/(EuclideanDistance[pt1, pt2]^2)];
inetforce[pt_, npt_] :=
Sum[force[pt, iposition[[i]]], {i, 1, npt}](*initial netforce*);
netforce[pt_, npt_] := Sum[force[pt, position[[i]]], {i, 1, npt}];
potential[position_, npt_] :=
Sum[Sum[1/EuclideanDistance[position[[i]], position[[j]]], {j,
i + 1, npt}], {i, 1, npt - 1}];
coord[\[Theta]_, \[Phi]_] := {Cos[\[Theta]] Sin[\[Phi]],
Sin[\[Theta]] Sin[\[Phi]], Cos[\[Phi]]};
newposition[p_, npt_] := Normalize[p + netforce[p, npt]];
\[Theta] = RandomReal[{0, 2 \[Pi]}, 470];
\[Phi] = RandomReal[{0, \[Pi]}, 470];
iposition =
Table[coord[\[Theta][[i]], \[Phi][[i]]], {i,
470}](*initial position*);
position = Table[0, 470];
Table[
Table[position[[i]] =
Normalize[iposition[[i]] + inetforce[iposition[[i]], n]], {i, n}];
Do[Table[position[[i]] = newposition[position[[i]], n] , {i, n}] ,
50 + 5 n];
Round[potential[position, n], .01]
, {n, 2, 20}]


In your code, you use Do loop to call the inside table 50+5n times. It seems like you want the potential value to converge. Based on this, the following is the tweak to check the potential value till it converges to 2 decimal places.

Firstly, define a function posR to call the table inside Do loop

posR[n_] :=
Table[position[[i]] = newposition[position[[i]], n], {i, n}];


Then, use a while loop to check whether the required precision is obtained or not.

Table[Table[
position[[i]] =
Normalize[iposition[[i]] + inetforce[iposition[[i]], n]], {i, n}];

While[Round[potential[posR[n], n], 0.001] =!=
Round[potential[posR[n], n], 0.001],
potential[posn[n], n];];

Round[potential[position, n], 0.01], {n, 2, 30}]//Timing
(*{2.23081, {0.5, 1.73, 3.67, 6.47, 9.99, 14.46, 19.68, 25.76, 32.72,
40.6, 49.17, 58.85, 69.31, 80.67, 92.91, 106.05, 120.09, 135.1,
150.88, 167.65, 185.3, 203.95, 223.35, 243.82, 265.16, 287.31,
310.5, 334.64, 359.61}}*)


The code you posted takes about 9.45 s to get the same result.

Although, there is an improvement by 4 times, it does not look optimised for higher values of n.

One suggestion would be to memoize some functions wherever applicable. Check this.

Hope this helps.

• Thanks! By using your code, the speed is definitely improved a lot. However, the values are still not accurate enough, I have to change the rounding precision to .00001 in While so that the final rounding is good enough. About using memorization, I think it can be used in force`, but it seems to slow down the calculation. – Taptic Feb 23 '17 at 23:29