# Distributions of Electrons on a Sphere Problem
# 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.
# See http://www-unix.mcs.anl.gov/~more/cops/.
gofn = function(dat, n)
{
x = dat[1:n]
y = dat[(n+1):(2*n)]
z = dat[(2*n+1):(3*n)]
ii = matrix(1:n, ncol = n, nrow = n, byrow = TRUE)
jj = matrix(1:n, ncol = n, nrow = n)
ij = which(ii<jj, arr.ind = TRUE)
i = ij[,1]
j = ij[,2]
# Coulomb potential
potential = sum(1.0/sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2))
potential
}
goeqfn = function(dat, n)
{
x = dat[1:n]
y = dat[(n+1):(2*n)]
z = dat[(2*n+1):(3*n)]
apply(cbind(x^2, y^2, z^2), 1, "sum")
}
n = 25
LB = rep(-1, 3*n)
UB = rep(1, 3*n)
eqB = rep(1, n)
ans = gosolnp(pars = NULL, fixed = NULL, fun = gofn, eqfun = goeqfn, eqB = eqB, LB = LB, UB = UB,
control = list(outer.iter = 100, trace = 1), distr = rep(1, length(LB)), distr.opt = list(),
n.restarts = 2, n.sim = 20000, rseed = 443, n = 25)
# should get a function value around 243.813
Run the code above in your browser using DataLab