Algorithm to place N points on a sphere with equal distances to all their closest neighbours. Uses the golden spiral algorithm to place the initial points followed by energy minimization using the steepest descent method to optimize the distances.

Steepest Descent

#!/usr/bin/env python

"""
Generates a set of points on a sphere and evenly distributes them
by minimizing the potential energy according to a 'Coloumb' potential.

Joao Rodrigues, 2013
"""

import numpy
from scipy.spatial import distance

def generate_points(N):
    """
    Generates a distribution of N points on a sphere based
    on the Golden Spiral Algorithm.
    """
    
    inc = numpy.pi * (3 - numpy.sqrt(5))
    off = 2.0 / N
    pts = [(numpy.cos(k * inc)*numpy.sqrt(1.0 - (k * off - 1 + (off / 2.0))**2), (k * off - 1 + (off / 2.0)), numpy.sin(k * inc)*numpy.sqrt(1.0 - (k * off - 1 + (off / 2.0))**2)) for k in xrange(N)]
    return numpy.array(pts, dtype=numpy.float64)

def calc_coulomb(coords):
    """
    Ecoulomb = E(1/r)
    """
    distances = distance.pdist(coords, 'euclidean')
    return numpy.power(distances, -1).sum()
    
def calc_forces(system):
    """
    Calculates forces acting on all particles.
    """
    
    forces = [ [0,0,0] for particle in system ]
    nparticles = len(system)
    
    for i in xrange(nparticles):
        particle_i = system[i]
        for j in xrange(i+1, nparticles):
            particle_j = system[j]
            r = particle_i - particle_j
            f = numpy.linalg.norm(r)**-3
            forces[i][0] += r[0]*f
            forces[i][1] += r[1]*f
            forces[i][2] += r[2]*f

            forces[j][0] -= r[0]*f
            forces[j][1] -= r[1]*f
            forces[j][2] -= r[2]*f
            
    return numpy.array(forces, dtype=numpy.float64)

def steepest_descent(coordinates, max_steps, step_size, min_step_size, min_energy_delta, verbose=True):
    """
    Steepest Descent minimization algorithm

    http://mathworld.wolfram.com/MethodofSteepestDescent.html
    """

    nparticles = len(coordinates)
    for k in xrange(options.nsteps):  
        # Backup coordinates
        old_coordinates = coordinates.copy()
        # 3.1 Calculate energy
        cur_energy = calc_coulomb(old_coordinates)

        # 3.2 Calculate forces
        forces = calc_forces(coordinates)
        force_normal = numpy.linalg.norm(forces)

        # 3.3 Move particles
        if force_normal > 0:
            for p in xrange(nparticles):
                directions = (forces[p]/force_normal)*step_size
                coordinates[p] += directions 
                coordinates[p] /= numpy.linalg.norm(coordinates[p])

        # 3.4 Check new energy and accept move if favourable + accelerate in that direction
        new_energy = calc_coulomb(coordinates)
        if new_energy >= cur_energy:
            step_size /= 2
            if step_size < min_step_size:
                print "-- Step Size Converged (%1.5f < %1.5f)" %(step_size, min_step_size)
                break
            coordinates = old_coordinates
        else:
            if verbose:
                print '-- Step %i: Energy %12.8f' %(k+1, cur_energy)
            delta_e = abs(new_energy - cur_energy)
            if delta_e < min_energy_delta:
                print "-- Energy Converged (%12.12f < %12.12f)" %(delta_e, min_energy_delta)
                break
            step_size *= 2
            cur_energy = new_energy
    
    return (coordinates, new_energy)

def plot_coordinates(initial, final):
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    x0, y0, z0 = zip(*initial)
    ax.scatter(x0, y0, z0, zdir='z', s=20, c='b')
    xmin, ymin, zmin = zip(*final)
    ax.scatter(xmin, ymin, zmin, zdir='z', s=20, c='r')
    u, v = numpy.mgrid[0:2*numpy.pi:20j, 0:numpy.pi:10j]
    x=numpy.cos(u)*numpy.sin(v)
    y=numpy.sin(u)*numpy.sin(v)
    z=numpy.cos(v)
    ax.plot_wireframe(x, y, z, color="gray")
    plt.show()
    
if __name__ == "__main__":
    
    import sys

    from optparse import OptionParser
    parser = OptionParser(usage="%prog <number of points> [options] [--help]")
    parser.add_option('-d', action="store", dest="step_size", type="float", 
                        default=0.01, help="Step size for the minimizer")
    parser.add_option('-n', action="store", dest="nsteps", type="int",
                        default=1000, help="Number of steps to minimize")
    parser.add_option('-c', action="store", dest="n_convergence", type="float",
                        default=(10**-5), help="Minimum Step Size (Convergence Criteria)")
    parser.add_option('-e', action="store", dest="e_convergence", type="float",
                        default=(10**-10), help="Minimum Energy Change (Convergence Criteria)")
    parser.add_option('-q', action="store_true", dest="quiet",
                        help="Supresses output from minimization algorithm")
    parser.add_option('-p', action="store_true", dest="plot",
                        help="Generates a plot of the initial and final coordinates [req. matplotlib]")
                        
    options, args = parser.parse_args()

    # 0. Input
    try:
        npoints = int(args[0])
    except (IndexError, ValueError):
        parser.print_usage()
        sys.exit(1)
    
    # 1. Generate Spiral Points
    initial_coordinates = generate_points(npoints)
    
    # 2. Minimization (steepest descent)
    optimized_coordinates, final_energy = steepest_descent(initial_coordinates.copy(), 
                                                            options.nsteps, 
                                                            options.step_size, 
                                                            options.n_convergence, 
                                                            options.e_convergence,
                                                            verbose=(not options.quiet))

    print '-- Final Energy: %12.8f' %(final_energy)

    # Plot initial vs final positions on the sphere
    if options.plot:
        plot_coordinates(initial_coordinates, optimized_coordinates)

    sys.exit(0)