# 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.