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.