Return average diameter for a molecule.

def find_average_diameter(elements, coordinates, adjust=1, increment=0.1,
                          processes=None, **kwargs):
    """Return average diameter for a molecule."""
    # Copy the coordinates as will perform many opertaions on them
    coordinates = deepcopy(coordinates)
    # Center of our cartesian system is always at origin
    origin = np.array([0, 0, 0])
    # Initial center of mass to reverse translation at the end
    initial_com = center_of_mass(elements, coordinates)
    # We just shift the cage to the origin.
    coordinates = shift_com(elements, coordinates)
    # We create an array of vdw radii of elements.
    elements_vdw = np.array([[atomic_vdw_radius[x.upper()]] for x in elements])
    # We calculate maximum diameter of a molecule to determine the radius
    # of a sampling sphere neccessary to enclose the whole molecule.
    shpere_radius = max_dim(elements, coordinates)[2]
    sphere_surface_area = 4 * np.pi * shpere_radius**2
    # Here we determine the number of sampling points necessary for a fine
    # sampling. Smaller molecules require more finner density of sampling
    # points on the sampling sphere's surface, whereas largen require less.
    # This formula was created so that larger molecule do not take much longer
    # to analyse, as number_sampling_points*length_of_sampling_vectors
    # results in quadratic increase of sampling time. The 250 factor was
    # specificly determined to produce close to 1 sampling point /Angstrom^2
    # for a sphere of radius ~ 24 Angstrom. We can adjust how fine is the
    # sampling by changing the adjust factor.
    number_of_points = int(np.log10(sphere_surface_area) * 250 * adjust)
    # Here I use code by Alexandre Devert for spreading points on a sphere:
    # http://blog.marmakoide.org/?p=1
    golden_angle = np.pi * (3 - np.sqrt(5))
    theta = golden_angle * np.arange(number_of_points)
    z = np.linspace(1 - 1.0 / number_of_points, 1.0 / number_of_points - 1.0,
                    number_of_points)
    radius = np.sqrt(1 - z * z)
    points = np.zeros((number_of_points, 3))
    points[:, 0] = radius * np.cos(theta) * shpere_radius
    points[:, 1] = radius * np.sin(theta) * shpere_radius
    points[:, 2] = z * shpere_radius
    # Here we analyse the vectors and retain the ones that create the molecule
    # outline.
    if processes:
        pool = Pool(processes=processes)
        parallel = [
            pool.apply_async(
                vector_analysis_reversed,
                args=(
                    point, coordinates, elements_vdw)
                ) for point in points
        ]
        results = [p.get() for p in parallel if p.get() is not None]
        pool.terminate()
    else:
        results = [
            vector_analysis_reversed(
                point, coordinates, elements_vdw)
            for point in points
        ]
    results_cleaned = [x[0] for x in results if x is not None]
    return np.mean(results_cleaned)*2