Python实现高密双尺寸球体堆积的优化方案问询

问题描述 #

我是Python编程新手，需要构建双尺寸球体的高密堆积结构，核心需求如下：

• 先让AM相（平均半径5）达到0.6的体积分数，再填充SE相（平均半径3）至0.25的体积分数，体积分数定义为对应相球体总体积与结构总体积的比值。

初始实现方案 #

采用随机生成球心的方法：生成坐标后检查重叠（允许一定程度重叠），满足条件则加入坐标列表；若球体跨边界，仅计算结构内部体积以统计体积分数，循环执行直至达到目标体积分数。

实现代码 #

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)

代码说明 #

• am_std_radius与se_std_radius用于实现尺寸正态分布，当前暂使用统一尺寸；
• check_overlap2函数用于控制球体的允许重叠限度；
• boundary_cross函数用于准确计算跨边界球体的内部体积，以正确统计体积分数。

当前困境 #

随机生成方法存在固有缺陷：当总体积分数达到约50%后，新粒子的生成耗时极长。希望更换为从底层构建高密结构的方案，了解到Drop and Roll等优化算法，但自行编码耗时久，需要实现思路或可用的GitHub仓库资源。

解决方案 #

实现思路（Drop and Roll算法） #

1. 初始密堆布局：先为AM相生成有序密堆积结构（比如面心立方FCC或六方密堆积HCP），这类结构理论最大体积分数约0.74，可轻松满足0.6的需求。
2. Drop阶段：对SE相小球，先随机生成候选球心，让小球沿重力方向（或指定方向）下落，直到碰到已有的AM相球体或结构边界。
3. Roll阶段：让下落到位的小球沿接触球体表面滚动，找到能量最低的稳定位置（尽可能贴合周围球体、减少空隙），同时保证不与其他球体过度重叠。
4. 体积分数校验：重复Drop和Roll步骤，定期统计已填充体积，直到SE相体积分数达到0.25。

可用的GitHub仓库资源 #

• Packmol封装库：Packmol是专业的粒子堆积工具，支持多尺寸球体密堆积，有Python封装库可直接调用，能快速生成符合体积分数要求的结构，无需自行实现复杂算法。
• pyscal：专注于粒子结构分析与生成的库，支持生成有序/无序密堆积结构，内置多种多尺寸粒子填充算法，文档清晰易上手。
• spheres：专门用于球体堆积结构生成的Python库，支持随机密堆积、有序密堆积及Drop and Roll优化算法，适合新手快速实现需求。

内容的提问来源于stack exchange，提问作者Sudarshan Wadajkar