zeogen_v8_refactored.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
RECIO-ZeoGen v8.0 (Refactored Object-Oriented Edition)
重构面向对象版本
[设计原则]
1. 参数与计算分离
2. 模块化设计
3. 可配置化接口
4. 易于扩展和维护
"""

import numpy as np
import random
from scipy.spatial import cKDTree
from collections import defaultdict
from abc import ABC, abstractmethod
from dataclasses import dataclass, field
from typing import List, Tuple, Dict, Optional, Any
from ase import Atoms
from ase.calculators.calculator import Calculator, all_changes
from ase.optimize import FIRE
from ase.io import write


# ==========================================
# 1. 配置参数模块
# ==========================================

@dataclass
class StructureParameters:
    """结构生成参数配置"""
    # 晶胞参数
    base_length_range: Tuple[float, float] = (10.0, 16.0)
    large_cell_probability: float = 0.1
    large_length_range: Tuple[float, float] = (20.0, 35.0)
    anisotropy_range: Tuple[float, float] = (0.7, 1.3)
    
    # 密度参数
    density_range: Tuple[float, float] = (15.0, 100.0)
    min_atoms: int = 4
    
    # 距离约束
    min_si_si_distance: float = 2.2
    sampling_max_attempts: int = 10000
    sampling_success_threshold: float = 0.5


@dataclass
class TopologyParameters:
    """拓扑参数配置"""
    # 傅里叶基函数配置
    n_basis_functions: int = 6
    basis_names: List[str] = field(default_factory=lambda: ['P', 'D', 'G', 'I-WP', 'Neo', '2X'])
    
    # 权重采样配置
    weight_sampling_mode: str = 'hybrid'  # 'pure', 'mixed', 'hybrid'
    pure_mode_probability: float = 0.3
    weight_alpha: List[float] = field(default_factory=lambda: [1.0, 1.0, 1.0, 1.0, 0.5, 0.5])
    
    # 阈值和各向异性
    threshold_range: Tuple[float, float] = (-0.8, 0.8)
    anisotropy_range: Tuple[float, float] = (0.8, 1.3)


@dataclass
class ForceFieldParameters:
    """力场参数配置"""
    # 键合参数
    si_o_bond_length: float = 1.61
    bond_force_constant: float = 15.0
    
    # 排斥参数
    repulsion_force_constant: float = 15.0
    max_cutoff_radius: float = 3.0
    
    # 物种特异性截断距离
    si_si_cutoff: float = 2.8
    o_o_cutoff: float = 2.2
    si_o_cutoff: float = 1.5
    
    # 力限制
    max_force_magnitude: float = 20.0


@dataclass
class OptimizationParameters:
    """优化参数配置"""
    optimizer_type: str = 'FIRE'
    max_force: float = 0.5
    max_steps: int = 120
    convergence_threshold: float = 0.5
    
    # 图构建参数
    initial_coordination: int = 3
    target_coordination: float = 4.0
    max_coordination: int = 6
    graph_balance_iterations: int = 5000
    neighbor_search_k: int = 24


@dataclass
class ZeoliteGenerationConfig:
    """总配置类"""
    structure: StructureParameters = field(default_factory=StructureParameters)
    topology: TopologyParameters = field(default_factory=TopologyParameters)
    force_field: ForceFieldParameters = field(default_factory=ForceFieldParameters)
    optimization: OptimizationParameters = field(default_factory=OptimizationParameters)
    
    # 输出配置
    output_format: str = 'cif'
    output_prefix: str = 'zeo_v8'
    
    def validate(self) -> bool:
        """验证参数合理性"""
        # 基础参数验证
        if self.structure.min_si_si_distance < 2.0:
            print("警告：Si-Si最小距离过小，可能导致结构不稳定")
        
        if self.force_field.si_o_bond_length < 1.5 or self.force_field.si_o_bond_length > 1.7:
            print("警告：Si-O键长偏离典型值")
        
        return True


# ==========================================
# 2. 力场计算模块
# ==========================================

class ForceFieldCalculator(ABC):
    """力场计算器抽象基类"""
    
    @abstractmethod
    def calculate_energy_and_forces(self, atoms: Atoms) -> Tuple[float, np.ndarray]:
        """计算能量和力"""
        pass


class ZeoliteDualRepulsionCalculator(Calculator, ForceFieldCalculator):
    """沸石双重排斥力场计算器"""
    
    def __init__(self, config: ForceFieldParameters, bonds_list: List[List[int]]):
        super().__init__()
        self.config = config
        self.bonds = np.array(bonds_list)
        self.bond_set = self._create_bond_set()
    
    def _create_bond_set(self) -> set:
        """创建键合对集合"""
        bond_set = set()
        if len(self.bonds) > 0:
            for b in self.bonds:
                bond_set.add(tuple(sorted((b[0], b[1]))))
        return bond_set
    
    def calculate(self, atoms=None, properties=['energy'], system_changes=all_changes):
        """ASE计算器接口"""
        super().calculate(atoms, properties, system_changes)
        energy, forces = self.calculate_energy_and_forces(atoms)
        self.results['energy'] = energy
        self.results['forces'] = forces
    
    def calculate_energy_and_forces(self, atoms: Atoms) -> Tuple[float, np.ndarray]:
        """计算能量和力的核心方法"""
        pos = atoms.get_positions()
        symbols = atoms.get_chemical_symbols()
        cell = atoms.get_cell()
        inv_cell = np.linalg.inv(cell)
        
        energy = 0.0
        forces = np.zeros_like(pos)
        
        # 键合能量
        energy += self._calculate_bond_energy(pos, forces, inv_cell, cell)
        
        # 排斥能量
        energy += self._calculate_repulsion_energy(pos, forces, symbols, inv_cell, cell)
        
        # 力限制
        self._apply_force_cap(forces)
        
        return energy, forces
    
    def _calculate_bond_energy(self, pos: np.ndarray, forces: np.ndarray, 
                             inv_cell: np.ndarray, cell: np.ndarray) -> float:
        """计算键合能量"""
        if len(self.bonds) == 0:
            return 0.0
        
        p1 = pos[self.bonds[:, 0]]
        p2 = pos[self.bonds[:, 1]]
        diff = p2 - p1
        d_frac = np.dot(diff, inv_cell)
        d_frac -= np.round(d_frac)
        diff_cart = np.dot(d_frac, cell)
        dists = np.linalg.norm(diff_cart, axis=1) + 1e-8
        
        delta = dists - self.config.si_o_bond_length
        energy = 0.5 * self.config.bond_force_constant * np.sum(delta**2)
        
        f_mag = -self.config.bond_force_constant * delta / dists
        f_vec = diff_cart * f_mag[:, np.newaxis]
        
        np.add.at(forces, self.bonds[:, 0], -f_vec)
        np.add.at(forces, self.bonds[:, 1], f_vec)
        
        return energy
    
    def _calculate_repulsion_energy(self, pos: np.ndarray, forces: np.ndarray,
                                  symbols: List[str], inv_cell: np.ndarray, 
                                  cell: np.ndarray) -> float:
        """计算排斥能量"""
        tree = cKDTree(pos)
        pairs = tree.query_pairs(r=self.config.max_cutoff_radius)
        
        energy = 0.0
        
        for i, j in pairs:
            if tuple(sorted((i, j))) in self.bond_set:
                continue
            
            s1, s2 = symbols[i], symbols[j]
            pair_type = tuple(sorted((s1, s2)))
            
            current_cutoff = self._get_cutoff_distance(pair_type)
            if current_cutoff == 0.0:
                continue
            
            # 计算距离
            d_vec = pos[j] - pos[i]
            d_frac = np.dot(d_vec, inv_cell)
            d_frac -= np.round(d_frac)
            d_vec = np.dot(d_frac, cell)
            r = np.linalg.norm(d_vec)
            if r < 0.05:
                r = 0.05
            
            if r < current_cutoff:
                energy += self._calculate_repulsion_term(r, current_cutoff, forces, i, j, d_vec)
        
        return energy
    
    def _get_cutoff_distance(self, pair_type: Tuple[str, str]) -> float:
        """获取截断距离"""
        cutoff_map = {
            ('Si', 'Si'): self.config.si_si_cutoff,
            ('O', 'O'): self.config.o_o_cutoff,
            ('O', 'Si'): self.config.si_o_cutoff
        }
        return cutoff_map.get(pair_type, 0.0)
    
    def _calculate_repulsion_term(self, r: float, cutoff: float, 
                                 forces: np.ndarray, i: int, j: int, 
                                 d_vec: np.ndarray) -> float:
        """计算排斥项"""
        scale = cutoff / r
        term = scale ** 4
        
        energy = self.config.repulsion_force_constant * term
        f_scalar = (4 * self.config.repulsion_force_constant * term / r**2)
        f_rep = f_scalar * d_vec
        
        forces[i] -= f_rep
        forces[j] += f_rep
        
        return energy
    
    def _apply_force_cap(self, forces: np.ndarray):
        """应用力限制"""
        f_norms = np.linalg.norm(forces, axis=1)
        large = f_norms > self.config.max_force_magnitude
        if np.any(large):
            forces[large] *= (self.config.max_force_magnitude / f_norms[large])[:, None]


# ==========================================
# 3. 拓扑生成模块
# ==========================================

class TopologyGenerator(ABC):
    """拓扑生成器抽象基类"""
    
    @abstractmethod
    def generate_surface_value(self, positions: np.ndarray, cell_dims: List[float],
                              coeffs: np.ndarray, threshold: float, 
                              anisotropy: np.ndarray) -> np.ndarray:
        """生成曲面值"""
        pass


class FourierSurfaceGenerator(TopologyGenerator):
    """傅里叶曲面生成器"""
    
    def __init__(self, config: TopologyParameters):
        self.config = config
    
    def generate_surface_value(self, positions: np.ndarray, cell_dims: List[float],
                              coeffs: np.ndarray, threshold: float, 
                              anisotropy: np.ndarray) -> np.ndarray:
        """生成傅里叶曲面值"""
        x, y, z = positions[:, 0], positions[:, 1], positions[:, 2]
        a, b, c = cell_dims
        
        # 基础频率
        kx, ky, kz = 2*np.pi/a, 2*np.pi/b, 2*np.pi/c
        sx, sy, sz = anisotropy
        X, Y, Z = kx*x*sx, ky*y*sy, kz*z*sz
        
        # 预计算三角函数
        cX, cY, cZ = np.cos(X), np.cos(Y), np.cos(Z)
        sX, sY, sZ = np.sin(X), np.sin(Y), np.sin(Z)
        c2X, c2Y, c2Z = np.cos(2*X), np.cos(2*Y), np.cos(2*Z)
        
        values = np.zeros(len(positions))
        
        # 基函数贡献
        if abs(coeffs[0]) > 1e-3:  # Primitive
            values += coeffs[0] * (cX + cY + cZ)
        
        if abs(coeffs[1]) > 1e-3:  # Diamond
            values += coeffs[1] * (cX*cY*cZ - sX*sY*sZ)
        
        if abs(coeffs[2]) > 1e-3:  # Gyroid
            values += coeffs[2] * (sX*cY + sY*cZ + sZ*cX)
        
        if abs(coeffs[3]) > 1e-3:  # I-WP
            term_iwp = 2*(cX*cY + cY*cZ + cZ*cX) - (c2X + c2Y + c2Z)
            values += coeffs[3] * term_iwp
        
        if abs(coeffs[4]) > 1e-3:  # Neovius
            term_neo = 3*(cX + cY + cZ) + 4*(cX*cY*cZ)
            values += coeffs[4] * term_neo
        
        if abs(coeffs[5]) > 1e-3:  # 2X
            values += coeffs[5] * (c2X + c2Y + c2Z)
        
        return values


class WeightSampler:
    """权重采样器"""
    
    def __init__(self, config: TopologyParameters):
        self.config = config
    
    def sample_weights(self, rng: np.random.RandomState) -> np.ndarray:
        """采样权重系数"""
        mode = self.config.weight_sampling_mode
        
        if mode == 'hybrid':
            current_mode = 'pure' if rng.rand() < self.config.pure_mode_probability else 'mixed'
        else:
            current_mode = mode
        
        if current_mode == 'pure':
            return self._sample_pure_weights(rng)
        elif current_mode == 'mixed':
            return self._sample_mixed_weights(rng)
        else:
            raise ValueError(f"Unknown sampling mode: {current_mode}")
    
    def _sample_pure_weights(self, rng: np.random.RandomState) -> np.ndarray:
        """采样纯权重"""
        n_bases = self.config.n_basis_functions
        w = np.zeros(n_bases)
        idx = rng.randint(0, n_bases)
        w[idx] = 1.0
        return w
    
    def _sample_mixed_weights(self, rng: np.random.RandomState) -> np.ndarray:
        """采样混合权重"""
        w = rng.dirichlet(self.config.weight_alpha)
        signs = rng.choice([-1, 1], size=self.config.n_basis_functions)
        return w * signs


# ==========================================
# 4. 结构构建模块
# ==========================================

class StructureBuilder:
    """结构构建器"""
    
    def __init__(self, config: OptimizationParameters):
        self.config = config
    
    def build_connectivity_graph(self, positions: np.ndarray, 
                                cell_matrix: np.ndarray) -> Tuple[Dict[int, set], Dict]:
        """构建连接图"""
        candidates, dist_map = self._get_neighbors(positions, cell_matrix, 
                                                  self.config.neighbor_search_k)
        
        # 目标边数
        target_edges = 2 * len(positions)
        adj = defaultdict(set)
        
        # 初始连接
        self._initial_connectivity(adj, candidates, dist_map, positions)
        
        # 动态平衡
        self._balance_graph(adj, candidates, dist_map, target_edges)
        
        return adj, dist_map
    
    def _initial_connectivity(self, adj: Dict[int, set], candidates: Dict, 
                           dist_map: Dict, positions: np.ndarray):
        """初始连接"""
        for i in range(len(positions)):
            pool = sorted(candidates[i], key=lambda x: dist_map.get((i, x), 999))
            for target in pool[:self.config.initial_coordination]:
                if target not in adj[i]:
                    adj[i].add(target)
                    adj[target].add(i)
    
    def _balance_graph(self, adj: Dict[int, set], candidates: Dict, 
                      dist_map: Dict, target_edges: int):
        """动态平衡图"""
        for _ in range(self.config.graph_balance_iterations):
            current_edges = sum([len(adj[i]) for i in adj]) // 2
            diff = current_edges - target_edges
            
            if diff == 0:
                break
            elif diff > 0:
                self._remove_edge(adj, dist_map)
            elif diff < 0:
                self._add_edge(adj, candidates, dist_map)
    
    def _remove_edge(self, adj: Dict[int, set], dist_map: Dict):
        """移除边"""
        high_deg = [i for i in adj if len(adj[i]) > self.config.max_coordination]
        if not high_deg:
            high_deg = list(adj.keys())
        
        node = random.choice(high_deg)
        if not adj[node]:
            return
        
        worst = max(adj[node], key=lambda x: dist_map.get(tuple(sorted((node, x))), 999))
        adj[node].remove(worst)
        adj[worst].remove(node)
    
    def _add_edge(self, adj: Dict[int, set], candidates: Dict, dist_map: Dict):
        """添加边"""
        low_deg = [i for i in range(len(adj)) if len(adj[i]) < self.config.target_coordination]
        if not low_deg:
            low_deg = list(range(len(adj)))
        
        node = random.choice(low_deg)
        pool = candidates[node]
        
        best = None
        min_d = 999.0
        for t in pool:
            if t not in adj[node]:
                d = dist_map.get((node, t), 999)
                if d < min_d:
                    min_d = d
                    best = t
        
        if best is not None:
            adj[node].add(best)
            adj[best].add(node)
    
    def _get_neighbors(self, pos: np.ndarray, cell_diag: np.ndarray, 
                      k: int) -> Tuple[Dict, Dict]:
        """获取邻居信息"""
        images = []
        for x in [-1, 0, 1]:
            for y in [-1, 0, 1]:
                for z in [-1, 0, 1]:
                    shift = np.dot([x, y, z], cell_diag)
                    images.append(pos + shift)
        
        super_pos = np.vstack(images)
        tree = cKDTree(super_pos)
        dists, idxs = tree.query(pos, k=k+1)
        
        n = len(pos)
        candidates = defaultdict(list)
        dist_map = {}
        inv = np.linalg.inv(cell_diag)
        
        for i in range(n):
            for d, super_idx in zip(dists[i], idxs[i]):
                orig_j = super_idx % n
                if i == orig_j:
                    continue
                
                vec = pos[orig_j] - pos[i]
                vf = np.dot(vec, inv)
                vf -= np.round(vf)
                real_d = np.linalg.norm(np.dot(vf, cell_diag))
                
                candidates[i].append(orig_j)
                dist_map[tuple(sorted((i, orig_j)))] = real_d
                dist_map[(i, orig_j)] = real_d
                dist_map[(orig_j, i)] = real_d
        
        return candidates, dist_map


# ==========================================
# 5. 优化模块
# ==========================================

class StructureOptimizer:
    """结构优化器"""
    
    def __init__(self, config: OptimizationParameters):
        self.config = config
    
    def optimize(self, atoms: Atoms, force_field: ForceFieldCalculator) -> bool:
        """优化结构"""
        try:
            atoms.calc = force_field
            dyn = FIRE(atoms, logfile=None)
            dyn.run(fmax=self.config.max_force, steps=self.config.max_steps)
            return True
        except Exception as e:
            print(f"优化失败: {e}")
            return False


# ==========================================
# 6. 主生成器类
# ==========================================

class ZeoliteGenerator:
    """沸石生成器主类"""
    
    def __init__(self, config: Optional[ZeoliteGenerationConfig] = None):
        self.config = config or ZeoliteGenerationConfig()
        self.config.validate()
        
        # 初始化组件
        self.topology_generator = FourierSurfaceGenerator(self.config.topology)
        self.weight_sampler = WeightSampler(self.config.topology)
        self.structure_builder = StructureBuilder(self.config.optimization)
        self.optimizer = StructureOptimizer(self.config.optimization)
    
    def generate(self, run_id: int) -> Tuple[Optional[Atoms], str]:
        """生成沸石结构"""
        rng = np.random.RandomState(run_id)
        
        try:
            # 1. 生成参数
            cell_params = self._sample_cell_parameters(rng)
            topology_params = self._sample_topology_parameters(rng)
            
            # 2. 采样Si位置
            si_positions = self._sample_si_positions(rng, cell_params, topology_params)
            if si_positions is None:
                return None, "Si位置采样失败"
            
            # 3. 构建连接图
            adj, dist_map = self.structure_builder.build_connectivity_graph(
                si_positions, cell_params['cell_matrix']
            )
            
            # 4. 组装原子结构
            atoms, bonds = self._assemble_atoms(si_positions, adj, cell_params)
            
            # 5. 几何优化
            force_field = ZeoliteDualRepulsionCalculator(
                self.config.force_field, bonds
            )
            success = self.optimizer.optimize(atoms, force_field)
            
            if not success:
                return None, "几何优化失败"
            
            # 6. 质量评估
            status = self._evaluate_structure(atoms, len(si_positions))
            
            atoms.wrap()
            return atoms, status
            
        except Exception as e:
            return None, f"生成过程出错: {str(e)}"
    
    def _sample_cell_parameters(self, rng: np.random.RandomState) -> Dict:
        """采样晶胞参数"""
        if rng.rand() < self.config.structure.large_cell_probability:
            L = rng.uniform(*self.config.structure.large_length_range)
        else:
            L = rng.uniform(*self.config.structure.base_length_range)
        
        cell_side = [L * rng.uniform(*self.config.structure.anisotropy_range) for _ in range(3)]
        cell = [*cell_side, 90, 90, 90]
        cell_matrix = np.diag(cell_side)
        
        return {
            'cell_side': cell_side,
            'cell': cell,
            'cell_matrix': cell_matrix,
            'density': rng.uniform(*self.config.structure.density_range)
        }
    
    def _sample_topology_parameters(self, rng: np.random.RandomState) -> Dict:
        """采样拓扑参数"""
        weights = self.weight_sampler.sample_weights(rng)
        threshold = rng.uniform(*self.config.topology.threshold_range)
        anisotropy = rng.uniform(*self.config.topology.anisotropy_range, 3)
        
        return {
            'weights': weights,
            'threshold': threshold,
            'anisotropy': anisotropy
        }
    
    def _sample_si_positions(self, rng: np.random.RandomState, 
                            cell_params: Dict, topology_params: Dict) -> Optional[np.ndarray]:
        """采样Si位置"""
        cell_side = cell_params['cell_side']
        cell_matrix = cell_params['cell_matrix']
        density = cell_params['density']
        
        n_target = int(np.prod(cell_side) / density)
        if n_target < self.config.structure.min_atoms:
            return None
        
        si_positions = []
        attempts = 0
        
        while len(si_positions) < n_target and attempts < self.config.structure.sampling_max_attempts:
            r = np.dot(rng.rand(3), cell_matrix)
            
            # 计算曲面值
            surface_val = self.topology_generator.generate_surface_value(
                r.reshape(1, 3), cell_side, 
                topology_params['weights'], 
                topology_params['threshold'],
                topology_params['anisotropy']
            )[0]
            
            # 接受概率
            if rng.rand() < np.exp(-(surface_val - topology_params['threshold'])**2 / 0.5):
                if not si_positions or self._min_distance(r, np.array(si_positions), cell_matrix) > self.config.structure.min_si_si_distance:
                    si_positions.append(r)
            
            attempts += 1
        
        if len(si_positions) < n_target * self.config.structure.sampling_success_threshold:
            return None
        
        return np.array(si_positions)
    
    def _min_distance(self, point: np.ndarray, others: np.ndarray, 
                     cell_matrix: np.ndarray) -> float:
        """计算最小距离"""
        inv = np.linalg.inv(cell_matrix)
        d = others - point
        df = np.dot(d, inv)
        df -= np.round(df)
        return np.min(np.linalg.norm(np.dot(df, cell_matrix), axis=1))
    
    def _assemble_atoms(self, si_positions: np.ndarray, adj: Dict[int, set], 
                       cell_params: Dict) -> Tuple[Atoms, List[List[int]]]:
        """组装原子结构"""
        n_si = len(si_positions)
        unique_edges = []
        
        for u in adj:
            for v in adj[u]:
                if u < v:
                    unique_edges.append((u, v))
        
        # 生成O位置
        o_positions = []
        bonds = []
        inv_cell = np.linalg.inv(cell_params['cell_matrix'])
        
        for u, v in unique_edges:
            p1 = si_positions[u]
            p2 = si_positions[v]
            diff = p2 - p1
            d_frac = np.dot(diff, inv_cell)
            d_frac -= np.round(d_frac)
            vec = np.dot(d_frac, cell_params['cell_matrix'])
            mid = p1 + 0.5 * vec
            o_positions.append(mid)
            
            o_idx = n_si + len(o_positions) - 1
            bonds.append([u, o_idx])
            bonds.append([v, o_idx])
        
        # 创建ASE原子对象
        atoms = Atoms(
            'Si' * n_si + 'O' * len(o_positions),
            positions=np.vstack([si_positions, np.array(o_positions)]),
            cell=cell_params['cell'],
            pbc=True
        )
        
        return atoms, bonds
    
    def _evaluate_structure(self, atoms: Atoms, n_si: int) -> str:
        """评估结构质量"""
        dists = atoms.get_all_distances(mic=True)
        si_cns = []
        
        for i in range(n_si):
            d = dists[i]
            cn = np.sum((d < 2.0) & (d > 0.1))
            si_cns.append(cn)
        
        avg_cn = np.mean(si_cns) if si_cns else 0
        formula = atoms.get_chemical_formula()
        
        status = "OK" if 3.5 < avg_cn < 4.5 else "DISTORTED"
        return f"{status} | {formula} | AvgCN={avg_cn:.2f}"


# ==========================================
# 7. 工厂类和便捷接口
# ==========================================

class ZeoliteGeneratorFactory:
    """沸石生成器工厂类"""
    
    @staticmethod
    def create_default_generator() -> ZeoliteGenerator:
        """创建默认生成器"""
        return ZeoliteGenerator()
    
    @staticmethod
    def create_large_pore_generator() -> ZeoliteGenerator:
        """创建大孔道生成器"""
        config = ZeoliteGenerationConfig()
        config.structure.density_range = (10.0, 30.0)
        config.topology.threshold_range = (-0.3, 0.2)
        return ZeoliteGenerator(config)
    
    @staticmethod
    def create_high_surface_area_generator() -> ZeoliteGenerator:
        """创建高比表面积生成器"""
        config = ZeoliteGenerationConfig()
        config.topology.weight_alpha = [0.5, 0.8, 0.5, 1.2, 0.3, 0.8]
        config.structure.density_range = (50.0, 120.0)
        return ZeoliteGenerator(config)
    
    @staticmethod
    def create_gyroid_dominant_generator() -> ZeoliteGenerator:
        """创建Gyroid主导生成器"""
        config = ZeoliteGenerationConfig()
        config.topology.weight_alpha = [0.1, 0.1, 2.0, 0.1, 0.1, 0.1]
        return ZeoliteGenerator(config)


# ==========================================
# 8. 主执行脚本
# ==========================================

def main():
    """主执行函数"""
    print(">>> RECIO-ZeoGen v8.0 (Refactored Object-Oriented Edition)")
    
    # 创建生成器
    generator = ZeoliteGeneratorFactory.create_default_generator()
    
    # 批量生成
    success_count = 0
    for i in range(20):
        atoms, message = generator.generate(i)
        
        if atoms:
            success_count += 1
            filename = f"{generator.config.output_prefix}_{success_count}.{generator.config.output_format}"
            write(filename, atoms)
            print(f"[{i+1}] {filename} -> {message}")
        else:
            print(f"[{i+1}] 生成失败: {message}")
    
    print(f"\n成功生成 {success_count}/20 个结构")


if __name__ == "__main__":
    main()