# Distribution of electrons on a sphere

TomSym implementation of GAMS Example (ELEC,SEQ=230)

Given n electrons, find the equilibrium state distribution (of minimal Coulomb potential) of the electrons positioned on a conducting sphere.

This model is from the COPS benchmarking suite.

The number of electrons can be specified using the parameter np.

Dolan, E D, and More, J J, Benchmarking Optimization Software with COPS. Tech. rep., Mathematics and Computer Science Division, 2000.

Morris, J R, Deaven, D M, and Ho, K M, Genetic Algorithm Energy Minimization for Point Charges on a Sphere. Phys. Rev. B. 53 (1996), R1740-R1743.

Saff, E B, and Kuijlaars, A, Distributing Many Points on the Sphere. Math. Intelligencer 19 (1997), 5-11.

```np = 25;

% Coordinates of the electrons
x = tom('x',np,1);
y = tom('y',np,1);
z = tom('z',np,1);

toms potential

% The following two methods of computing the potential energy are
% mathematically equivalent.
useSlowMethod = false;  % Set to "true" to try the slower method. (Tip: Use
% a small value for np when you try this.)

if useSlowMethod
% The slow way of computing the sum
% Using nested for loops, we create a symbolic expression containing
% (np-1)*(np-2)/2 terms. Each of these terms will be symbolically
% differentiated to compute the Jacobian and Hessian matrices for use
% by the solver.
% For np=10, the autogenerated Jacobian of this function contains more
% than 15 kilobytes of Matlab code (and the Hessian more than 30 kb).
% TomSym will need several minutes to generate this code, which is
% orders of magnitude more than the fraction of a second required to
% actually solve the numeric problem.
potential = 0;
for i=1:np
for j=(i+1):np
potential = potential + 1.0/sqrt((x(i)-x(j))^2 +...
(y(i)-y(j))^2 + (z(i)-z(j))^2);
end
end
else
% Using vectorized code saves a lot of time in handling the symbolic
% expressions, because the symbolic derivative is computed for the one
% vectorized expression, not on hundreds of individual terms.
% (The autogenerated Hessian of this function contains less than 1 kb
% of Matlab code, regardless of the value of np.)
[ii,jj] = meshgrid(1:np,1:np);
[i,j] = find(ii<jj);
potential = sum(1.0./sqrt((x(i)-x(j)).^2 + (y(i)-y(j)).^2 +...
(z(i)-z(j)).^2));
end

% Ball constraint
ball = (x.^2 + y.^2 + z.^2 == 1);

% Initial guess:
theta = 2*pi*rand(np,1);
phi   = pi*rand(np,1);
x0 = {x == cos(theta).*sin(phi)
y == sin(theta).*sin(phi)
z == cos(phi)};

options = struct;
options.solver = 'snopt';
solution = ezsolve(potential,ball,x0,options);
```
```Problem type appears to be: con
===== * * * =================================================================== * * *
TOMLAB - Tomlab Optimization Inc. Development license  999001. Valid to 2010-02-05
=====================================================================================
Problem: ---  1: Problem 1                      f_k     243.812760298469870000
sum(|constr|)      0.000000000064759309
f(x_k) + sum(|constr|)    243.812760298534640000
f(x_0)    326.435082121151590000

Solver: snopt.  EXIT=0.  INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv  174 GradEv  172 ConstrEv  172 ConJacEv  172 Iter  171 MinorIter  224
CPU time: 0.484375 sec. Elapsed time: 0.484000 sec.
```