Subject: Re: CCL:charges on sphere, summary

From: elewars <[email-masked]>
Date: Wed, 10 Oct 2001 21:18:23 -0400

Hello,

Thank you for your very interesting contribution to this discussion; and your idea of _simulating_ the situation is
novel.
I am passing on your response to the CCL net and to my colleague.

E. Lewars
=====

Ivan Tubert wrote:

> Hi,
>
> I have a new response about the charge problem.
>
> Qualitatively, I think that they should all, or maybe all but one, be on
> the surface.
>
> * We know that for small n, such as n=2, they are all on the surface.
>
> * In the limit n->infinity, if we assume that we have the charges on the
> surface and add a new one, we have that the potential energy inside the
> sphere is constant (it is the classical example of a charged sphere).
>
> * Therefore, when adding a new charge it would make no difference whether
> we add it to the surface or inside the sphere. Let's assume we add it
> to the center of the sphere, just for the sake of argument. However, we
> cannot have two charges inside the sphere because they would repel each
> other just as in the case n=2.
>
> So, for very large n, we would have n-1 charges on the surface, and one
> "free"  charge.
>
> This does not look like a formal proof, of course. So I decided to make an
> "experiment" with a simple computer simulation. I wrote a program to
> optimize the energy of the system by simulated annealing, and the results
> up to n=100 always end up with all of the charges on the surface. However,
> for large n, the last few charges (and especially the last one) take a
> while to migrate to the surface, seemingly in agreement with the concept
> of "one free charge".
>
> The program is written in C, is not user-friendly, and surely it is very
> inefficient (I would not recommend trying to optimize n > 100). You
> have to edit the source and compile it to change the parameters, but I
> include the source as an attachment in case you are interested.
>
> Regards,
>   Ivan Tubert
>
>   ------------------------------------------------------------------------
>
> #include <stdlib.h>
>
> #include <stdio.h>
>
> #include <string.h>
>
> #include <math.h>
>
> #include <assert.h>
>
> #include <time.h>
>
> typedef struct {
>
>         double x, y, z;
>
>         double q;
>
> } PARTICLE;
>
> typedef struct {
>
>         double r, theta, phi;
>
> } POLAR;
>
> #define PI 3.1415926535
>
> /************* PARAMETERS ***********/
>
> double r = 1.0;                         // radius of the sphere
>
> double k = 1.0e-4;                      // "Boltzmann constant"
>
> double stepsize = 0.005;        // size of random move
>
> double T0 = 1000.0;                     // initial temperature
>
> double DT = 50.0;                       // temperature step
>
> int NSTEPS = 5000;                      // number of moves for each temperature
>
> int NPART = 48;                         // Number of particles
>
> /************************************/
>
> double get_r(PARTICLE p)
>
> {
>
>         return(sqrt(p.x*p.x+p.y*p.y+p.z*p.z));
>
> }
>
> PARTICLE *init(int n, double q)
>
> {
>
>         int i;
>
>         PARTICLE *p;
>
>         p = malloc(sizeof(PARTICLE)*n);
>
>         assert (p);
>
>         for(i = 0; i < n; ++i){
>
>                 p[i].q = q;
>
>                 do {
>
>                         p[i].x = (1.0*rand()/RAND_MAX - 0.5) * 2 * r;
>
>                         p[i].y = (1.0*rand()/RAND_MAX - 0.5) * 2 * r;
>
>                         p[i].z = (1.0*rand()/RAND_MAX - 0.5) * 2 * r;
>
>                 } while (get_r(p[i])>r);
>
>         }
>
>
>
>         return p;
>
> }
>
> double get_dist(PARTICLE p1,PARTICLE p2)
>
> {
>
>         return( sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z)) );
>
> }
>
> double get_E(PARTICLE *p)
>
> {
>
>         int i,j;
>
>         double E = 0.0;
>
>         for(i = 0; i < NPART-1; ++i){
>
>                 for(j = i+1; j < NPART; ++j){
>
>                         E += p[i].q*p[j].q/get_dist(p[i], p[j]);
>
>                 }
>
>         }
>
>         return E;
>
> }
>
> POLAR get_polar(PARTICLE p)
>
> {
>
>         POLAR pol;
>
>         pol.r = get_r(p);
>
>         pol.theta = acos(p.z/pol.r);
>
>         pol.phi = atan(p.y/p.x);
>
>         if(p.y < 0)
>
>                 pol.phi += PI/2;
>
>         return pol;
>
> }
>
> void get_stat_r(PARTICLE *p, double *mean, double *sd)
>
> {
>
>         int i,j;
>
>         double r = 0.0, r2 = 0.0, d;
>
>
>
>         for(i = 0; i < NPART-1; ++i){
>
>                 for(j = i+1; j < NPART; ++j){
>
>                         d = get_dist(p[i], p[j]);
>
>                         r += d;
>
>                         r2 += d*d;
>
>                 }
>
>         }
>
>         *mean = 2*r/NPART/(NPART-1);
>
>         *sd = sqrt(2*r2/NPART/(NPART-1) - *mean * *mean);
>
> }
>
> void print(PARTICLE *p)
>
> {
>
>         int i;
>
>         POLAR pol;
>
>         double mean, sd;
>
>         printf("     x        y      z   ->      r     theta     phi\n");
>
>         for(i = 0; i < NPART; ++i){
>
>                 pol = get_polar(p[i]);
>
>                 printf("%8.2f%8.2f%8.2f -> %8.2f%8.2f%8.2f\n", p[i].x, p[i].y, p[i].z, pol.r, pol.theta, pol.phi);
>
>         }
>
>
>
>         get_stat_r(p, &mean, &sd);
>
>         printf("    mean radius r(ij) = %5.2f +/- %5.2f\n", mean, sd);
>
> }
>
> int main()
>
> {
>
>         double T;
>
>         int i;
>
>         PARTICLE *p, newp, oldp;
>
>         double oldE, newE;
>
>         int ip;
>
>         double bz,rn;
>
>         int reject;
>
>
>
>         srand( (unsigned int)time( NULL ) );
>
>         T = T0;
>
>         p = init(NPART, 1.0);
>
>         newp.q = 1.0;
>
>         print(p);
>
>         oldE = get_E(p);
>
>         while(T >= 0.0) {
>
>                 reject = 0;
>
>                 printf("T = %5.f ... ", T);
>
>                 for(i = 0; i < NSTEPS; ++i){
>
>                         ip = rand()%NPART;
>
>                         do {  // move
>
>                                 newp.x = p[ip].x + (1.0*rand()/RAND_MAX - 0.5) * 2 * stepsize;
>
>                                 newp.y = p[ip].y + (1.0*rand()/RAND_MAX - 0.5) * 2 * stepsize;
>
>                                 newp.z = p[ip].z + (1.0*rand()/RAND_MAX - 0.5) * 2 * stepsize;
>
>                         } while (get_r(newp)>r);
>
>                         oldp = p[ip];
>
>                         p[ip] = newp;
>
>                         newE = get_E(p);
>
>                         if(newE > oldE){
>
>                                 bz = exp((oldE-newE)/k/T);
>
>                                 rn = 1.0*rand()/RAND_MAX;
>
>                                 if (rn > bz) { // reject
>
>                                         p[ip] = oldp;
>
>                                         newE = oldE;
>
>                                         ++reject;
>
>                                 }
>
>                         }
>
>
>
>                         oldE = newE;
>
>                 }
>
>                 T -= DT;
>
>                 printf("Rejected %i moves out of %i\n", reject, NSTEPS);
>
>                 print(p);
>
>         }
>
>
>
>         return(0);
>
> }