Return to Answer

pts = Apply[{2 \[Pi] #1, ArcCos[2 #2 - 1]} &, RandomReal[1, {10, 2}], 1];
Clear[energy];
Clear[a];
vars = Array[a, {Length[pts], 2}];
energy[p_] := 
  Module[{cart}, 
   cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]], 
       Cos[#[[1]]]} &, p];
   Total[Outer[Exp[-Norm[#1 - #2]] &, cart, cart, 1], 2]];
FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts}]]

pts = Apply[{2 \[Pi] #1, ArcCos[2 #2 - 1]} &, RandomReal[1, {100, 2}], 1];
Clear[energy];
Clear[a];
vars = Array[a, {Length[pts], 2}];
energy[p_] := 
  Module[{cart}, 
   cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]], 
       Cos[#[[1]]]} &, p];
   Total[Outer[Exp[-Sqrt[(#1 - #2).(#1 - #2)]] &, cart, cart, 1], 2]];
FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts}]]

pts = Apply[{ArcCos[2 #2 - 1], 2 \[Pi] #1} &, RandomReal[1, {10, 2}], 1];
Clear[energy];
Clear[a];
vars = Array[a, {Length[pts], 2}];
energy[p_] := 
  Module[{cart}, 
   cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]], 
       Cos[#[[1]]]} &, p];
   Total[Outer[Exp[-Norm[#1 - #2]] &, cart, cart, 1], 2]];
FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts}]]

pts = Apply[{ArcCos[2 #2 - 1], 2 \[Pi] #1} &, RandomReal[1, {100, 2}], 1];
Clear[energy];
Clear[a];
vars = Array[a, {Length[pts], 2}];
energy[p_] := 
  Module[{cart}, 
   cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]], 
       Cos[#[[1]]]} &, p];
   Total[Outer[Exp[-Sqrt[(#1 - #2).(#1 - #2)]] &, cart, cart, 1], 2]];
FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts}]]

Oh, and one more thing I changed (unrelated to the question), is to switch your definitions of polar and azimuthal angles in pts.

Edit 4

I just had another idea on how to improve the speed of my solution: the use of Norm might make it harder to estimate the Hessian for this function. And indeed, when I got rid of Norm there was a significant speed gain (note that the initial solution above is already faster than the _?NumericQ approach even when the latter is compiled while mine is not). I think this is worth adding here because Norm seems like a natural thing to use in pair potentials, even if the energy expression becomes more complicated than the one in this question.

So here is the new version, with Norm replaced by Sqrt[(#1 - #2).(#1 - #2)]. Observe that I have now put back the original particle number of 100 because on my laptop this takes less than 8 seconds to evaluate!

pts = Apply[{2 \[Pi] #1, ArcCos[2 #2 - 1]} &, RandomReal[1, {100, 2}], 1];
Clear[energy];
Clear[a];
vars = Array[a, {Length[pts], 2}];
energy[p_] := 
  Module[{cart}, 
   cart = Map[{Sin[#[[1]]]*Cos[#[[2]]], Sin[#[[1]]]*Sin[#[[2]]], 
       Cos[#[[1]]]} &, p];
   Total[Outer[Exp[-Sqrt[(#1 - #2).(#1 - #2)]] &, cart, cart, 1], 2]];
FindMinimum[energy[vars], Transpose[{Flatten@vars, Flatten@pts}]]

I get a significant speed-up with this for your example, but the performance depends on the random starting points (and on the choice of bracket width) so I can't say anything definitive. That seems like a topic for a different question.

I did get a significant speed-up with this for your example, but the performance depends on the random starting points (and on the choice of bracket width) so I can't say anything definitive. That seems like a topic for a different question.

Edit 3

Though I didn't look at the speed issue in detail, forcing FindMinimum to work with numerical derivatives may be the worst option here. That will happen if you define your function energy only for numerical arguments, as in

energy[p : {{_?NumericQ, _?NumericQ} ..}] := 

followed by either your own or my initial definition above. So although that's a common advice people give in these applications, it is not going to be the fastest approach here.