Creating a dense packing of spheres

Question

I am new to python coding and I am having a very hard time to create a dense packing of spheres.

Problem and my approach towards it:

I am creating a structure filled with spheres of two different sizes. Let's say they are named phase 1 with radius=r1 and phase 2 with radius=r2. I want phase 1 to reach a volume fraction of vf1 and then fill the structure with phase 2 spheres until they reach volume fraction vf2. Volume fraction is the ratio of total volume of spheres(respective phase) to the total volume of the structure.

** Flow -> ** So I generate a coordinate, check the overlap, if overlap check satisfied then add to our coordinates list, if it is crossing boundary then consider only the volume inside for volume fraction, repeat.

Initially I started with choosing random centers in my structure and generated particles until a given volume fraction is reached.

**Here is the code(My naming of the phases might be different but the idea is the same): **

import numpy as np
import time
from tqdm import tqdm

import numpy as np
from skimage import io, draw

# Define microstructure parameters
size = (100, 100, 100)  # Microstructure dimensions
am_mean_radius = 5  # Mean AM particle radius
deviation_fraction=0.1
am_std_radius = am_mean_radius*deviation_fraction  # Standard deviation of AM particle radius
se_mean_radius = 3  # Mean SE particle radius
se_std_radius = 0.5  # Standard deviation of SE particle radius
am_fraction = 0.6  # Desired volume fraction of the AM particles
se_fraction = 0.25  # Desired volume fraction of the SE particles
max_range = 2 * am_mean_radius
min_distance = 0
overlap_fraction=0.9 # This value multiplied to R1+R2 


# Initialize the particle coordinates
am_coordinates = []
se_coordinates = []

#Intialize the particle volume
am_vol=0 # Total AM particle volume
se_vol=0 # Total SE particle volume

def check_overlap2(coordinates, center, radius):    
    if not coordinates:
        return True
    
    # Extract the coordinates within the cube defined by the radius
    x_range = (center[0] - max_range, center[0] + max_range)
    y_range = (center[1] - max_range, center[1] + max_range)
    z_range = (center[2] - max_range, center[2] + max_range)
    
    filtered_coordinates = [(c[0], c[1]) for c in coordinates if
                        x_range[0] <= c[0][0] <= x_range[1] and
                        y_range[0] <= c[0][1] <= y_range[1] and
                        z_range[0] <= c[0][2] <= z_range[1]]

    for existing_center, existing_radius in filtered_coordinates:
        distance = np.linalg.norm(np.array(existing_center) - np.array(center))
        if distance < (overlap_fraction) * (existing_radius + radius):
            return False
            
    return True


# Check if a particle crosses the microstructure boundary.
def boundary_cross(center, radius, microstructure_size):
    """
    Check if a particle crosses the microstructure boundary.

    Parameters:
        center (numpy.ndarray): Center coordinates of the particle.
        radius (float): Radius of the particle.
        microstructure_size (numpy.ndarray): Size of the microstructure.

    Returns:
        bool: True if the particle crosses the boundary, False otherwise.
    """
    return np.any(center < radius) or np.any(center >= np.array(microstructure_size) - radius)



# Function to generate alternating AM and SE particles with normal distribution of sizes
def generate_structure(size, am_mean_radius, am_std_radius, se_mean_radius, se_std_radius, min_distance, am_fraction, se_fraction):
    am_coordinates = []
    se_coordinates = []
    coordinates = []
    am_vol = 0
    se_vol = 0
    component = "AM"  # Starting with AM component
    target_am_vol = am_fraction * np.prod(size)
    target_se_vol = se_fraction * np.prod(size)
    target_vol = target_am_vol + target_se_vol
    progress_bar = tqdm(total=target_vol, desc="Generating Particles")
    
    # Set the time limit in seconds
    time_limit = 7200  # 0.5 hour
    start_time = time.time()
    elapsed_time = 0
    
    while am_vol < target_am_vol or se_vol < target_se_vol:
        
        if elapsed_time >= time_limit:
            print("Time limit exceeded. Exiting the loop.")
            break
            
        if component == "AM":
            radius = am_mean_radius
            vol = am_vol
            target_component_vol = target_am_vol
        else:
            radius = se_mean_radius
            vol = se_vol
            target_component_vol = target_se_vol

        while True:
            center = np.random.rand(3) * np.array(size)

            if check_overlap2(coordinates, center, radius):
                if boundary_cross(center, radius, size):
                    d = np.min([np.abs(size - center), np.abs(center)])
                    intersected_height = radius - d
                    intersected_volume = (1/3) * np.pi * intersected_height**2 * (3 * radius - intersected_height)
                    complete_volume = (4/3) * np.pi * radius**3
                    inside_volume = complete_volume - intersected_volume
                    vol = vol + inside_volume
                else:
                    vol = vol + 4/3 * np.pi * radius**3

                coordinates.append((center, radius, component))
                progress_bar.update((am_vol+se_vol) - progress_bar.n)  # Update progress bar based on the added volume
                # Update the elapsed time
                elapsed_time = time.time() - start_time
                break

        if component == "AM":
            am_vol = vol
            if am_vol >= target_am_vol:
                component = "SE"  # Switch to SE component
        else:
            se_vol = vol
            if se_vol >= target_se_vol:
                component = "AM"  # Switch to AM component
        
        
        
    progress_bar.close()  # Close the progress bar

    for center, radius, particle_type in coordinates:
        if particle_type == 'AM':
            am_coordinates.append((center, radius))
        else:
            se_coordinates.append((center, radius))

    return am_coordinates, se_coordinates, am_vol, se_vol, coordinates


# Generate alternating AM and SE particles with normal distribution of sizes
am_coordinates, se_coordinates, am_vol, se_vol, coordinates = generate_structure(size, am_mean_radius, am_std_radius, se_mean_radius, se_std_radius, min_distance, am_fraction, se_fraction)

Note:

1. am_std_radius and se_std_radius are for applying normal distribution of sizes but for the start I am applying it for uniform sized (Apologies if I have introduced some terms abruptly).
2. check_overlap2 function is to allow overlap within a certain limit.
3. boundary_cross function helps to correctly calculate volume fraction as some spheres might cross the boundary and thus to calculate them correctly we use it here.

Now my problem: The structure can be created for a total volume fraction of around 50% but after that generating every new particle generated is taking eternity and it should so because I am just generating it randomly. I want to change the approach and want to create a dense structure from scratch. I had looked for some optimized algorithms on the internet such as drop and roll algorithm for spherical packing but it will take a lot of time for me to code the same. So it would be great if you could tell me how to go on about this problem or maybe if you have some github repositories that would help a lot too.