""" TLT-006: Blind N-Fold Assignment Validation ============================================= Tests whether the N-fold symmetry mapping used in TLT-004 can be independently validated by examining 2D projections of 3D crystal structures along all high-symmetry directions. AUDIT CONTEXT: Cross-AI audit (Entry 9) scored N-fold mapping 2/10 (Gemini), 6/10 (Grok). The mapping (FCC→3, HCP→3, BCC→4, diamond→3, SC→4) was called "arbitrary" and "ad hoc." This test validates it empirically rather than deriving it from first principles. METHOD: 1. For each standard crystal type (FCC, BCC, HCP, diamond, SC, rhombo, BCT): a. Generate a supercell of 3D lattice points b. Project onto planes perpendicular to all high-symmetry directions c. For each projection, compute: - 2D planar density (atoms per unit area) - Nearest-neighbor angular distribution - Rotational symmetry order (2, 3, 4, or 6-fold) 2. The "blind N" is the symmetry order of the DENSEST projection plane 3. Compare blind N against model assignment 4. Report concordance/discordance CRITERION: The densest crystallographic plane is an objective, physical property — it determines cleavage, growth, slip directions, and surface stability. Using it as the selection criterion is physically motivated, not arbitrary. For cubic systems: [100], [110], [111] projections tested For hexagonal systems: [0001], [10-10], [11-20] projections tested Date: 2026-03-13 Test ID: TLT-006 """ import numpy as np import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from pathlib import Path from collections import Counter OUTPUT_DIR = Path("/root/rhetroiluso/project_prometheus/time_ledger_theory/" "tlt results/unaudited/TLT-006_blind_N") OUTPUT_DIR.mkdir(parents=True, exist_ok=True) # The MODEL's assignments (what we're testing against) MODEL_MAPPING = { 'FCC': 3, 'HCP': 3, 'BCC': 4, 'diamond': 3, 'SC': 4, 'BCT': 4, 'rhombo': 3, } # ============================================================================ # CRYSTAL LATTICE GENERATION # ============================================================================ def generate_supercell(structure, n_cells=4): """Generate 3D lattice points for standard crystal structures. Uses normalized lattice parameter a=1.0. Returns Cartesian coordinates. """ a = 1.0 if structure == 'FCC': basis_frac = np.array([ [0.0, 0.0, 0.0], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5], ]) cell = np.eye(3) * a elif structure == 'BCC': basis_frac = np.array([ [0.0, 0.0, 0.0], [0.5, 0.5, 0.5], ]) cell = np.eye(3) * a elif structure == 'SC': basis_frac = np.array([[0.0, 0.0, 0.0]]) cell = np.eye(3) * a elif structure == 'diamond': basis_frac = np.array([ [0.0, 0.0, 0.0], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5], [0.25, 0.25, 0.25], [0.75, 0.75, 0.25], [0.75, 0.25, 0.75], [0.25, 0.75, 0.75], ]) cell = np.eye(3) * a elif structure == 'HCP': c = a * np.sqrt(8.0 / 3.0) # ideal c/a a1 = np.array([a, 0, 0]) a2 = np.array([-a / 2, a * np.sqrt(3) / 2, 0]) a3 = np.array([0, 0, c]) cell = np.array([a1, a2, a3]) basis_frac = np.array([ [0.0, 0.0, 0.0], [1./3., 2./3., 0.5], ]) elif structure == 'BCT': # Body-centered tetragonal (c/a = 1.2, between FCC and BCC) c = a * 1.2 basis_frac = np.array([ [0.0, 0.0, 0.0], [0.5, 0.5, 0.5], ]) cell = np.diag([a, a, c]) elif structure == 'rhombo': # Rhombohedral (alpha = 60 degrees, like Bi/Sb A7 structure) alpha = np.radians(60) cos_a = np.cos(alpha) sin_a = np.sin(alpha) # Rhombohedral cell vectors a1 = a * np.array([1, 0, 0]) a2 = a * np.array([cos_a, sin_a, 0]) a3_z = np.sqrt(1 - 2 * cos_a**2 + cos_a) if cos_a < 1 else 0 a3 = a * np.array([cos_a, (cos_a - cos_a**2) / sin_a if sin_a > 0 else 0, np.sqrt(1 - cos_a**2 - ((cos_a - cos_a**2) / sin_a)**2) if sin_a > 0 else 1]) cell = np.array([a1, a2, a3]) basis_frac = np.array([[0.0, 0.0, 0.0]]) else: raise ValueError(f"Unknown structure: {structure}") # Generate supercell points = [] n = n_cells for i in range(-n, n + 1): for j in range(-n, n + 1): for k in range(-n, n + 1): translation = i * cell[0] + j * cell[1] + k * cell[2] for bf in basis_frac: cart = bf[0] * cell[0] + bf[1] * cell[1] + bf[2] * cell[2] points.append(cart + translation) return np.array(points), cell def get_projection_directions(structure): """Return high-symmetry projection directions for each crystal system.""" if structure in ('FCC', 'BCC', 'SC', 'diamond', 'BCT'): return { '[100]': np.array([1, 0, 0], dtype=float), '[110]': np.array([1, 1, 0], dtype=float) / np.sqrt(2), '[111]': np.array([1, 1, 1], dtype=float) / np.sqrt(3), } elif structure == 'HCP': a = 1.0 c = a * np.sqrt(8.0 / 3.0) return { '[0001]': np.array([0, 0, 1], dtype=float), '[10-10]': np.array([1, 0, 0], dtype=float), '[11-20]': np.array([1, np.sqrt(3) / 2, 0], dtype=float), } elif structure == 'rhombo': return { '[111]': np.array([1, 1, 1], dtype=float) / np.sqrt(3), '[110]': np.array([1, 1, 0], dtype=float) / np.sqrt(2), '[100]': np.array([1, 0, 0], dtype=float), } else: return { '[100]': np.array([1, 0, 0], dtype=float), '[110]': np.array([1, 1, 0], dtype=float) / np.sqrt(2), '[111]': np.array([1, 1, 1], dtype=float) / np.sqrt(3), } # ============================================================================ # PROJECTION AND SYMMETRY ANALYSIS # ============================================================================ def project_onto_plane(points_3d, normal): """Project 3D points onto a 2D plane perpendicular to normal.""" normal = normal / np.linalg.norm(normal) # Find two orthogonal basis vectors in the plane if abs(normal[2]) < 0.9: v1 = np.cross(normal, np.array([0, 0, 1])) else: v1 = np.cross(normal, np.array([1, 0, 0])) v1 = v1 / np.linalg.norm(v1) v2 = np.cross(normal, v1) v2 = v2 / np.linalg.norm(v2) # Project coords_2d = np.column_stack([ np.dot(points_3d, v1), np.dot(points_3d, v2), ]) return coords_2d def deduplicate_2d(coords_2d, tol=0.01): """Remove duplicate 2D points (within tolerance).""" unique = [coords_2d[0]] for pt in coords_2d[1:]: dists = np.linalg.norm(np.array(unique) - pt, axis=1) if np.min(dists) > tol: unique.append(pt) return np.array(unique) def extract_single_plane(points_3d, normal, tol=0.05): """Extract atoms lying on a single crystallographic plane. Instead of projecting ALL atoms (which conflates layers), this finds the distinct (hkl) planes perpendicular to `normal` and returns atoms from the most populated plane near the origin. Returns: (plane_atoms_2d, d_spacing, n_layers) """ normal = normal / np.linalg.norm(normal) # Signed distance of each atom along the normal signed_dists = np.dot(points_3d, normal) # Cluster into layers by signed distance sorted_dists = np.sort(np.unique(np.round(signed_dists, decimals=3))) layers = [] used = np.zeros(len(signed_dists), dtype=bool) for d in sorted_dists: mask = (~used) & (np.abs(signed_dists - d) < tol) if np.sum(mask) > 0: layers.append((d, np.where(mask)[0])) used |= mask if len(layers) == 0: return None, 0, 0 # Find the d-spacing (distance between adjacent layers) layer_ds = sorted([l[0] for l in layers]) if len(layer_ds) > 1: spacings = np.diff(layer_ds) d_spacing = np.median(spacings[spacings > tol]) else: d_spacing = 0 # Pick the layer closest to origin with enough atoms best_layer = None best_dist = float('inf') for d, indices in layers: if len(indices) >= 3 and abs(d) < best_dist: best_dist = abs(d) best_layer = (d, indices) if best_layer is None: # Fallback: pick the most populated layer best_layer = max(layers, key=lambda l: len(l[1])) _, indices = best_layer plane_atoms_3d = points_3d[indices] # Project these atoms onto the 2D plane plane_atoms_2d = project_onto_plane(plane_atoms_3d, normal) plane_atoms_2d = deduplicate_2d(plane_atoms_2d, tol=0.03) return plane_atoms_2d, d_spacing, len(layers) def compute_planar_density(unique_2d, radius=3.0): """Compute planar density: points per unit area within radius.""" center = np.mean(unique_2d, axis=0) dists = np.linalg.norm(unique_2d - center, axis=1) n_in_circle = np.sum(dists <= radius) area = np.pi * radius**2 return n_in_circle / area def compute_angular_symmetry(unique_2d, k_neighbors=6, radius=2.5): """Compute the dominant rotational symmetry order from nearest-neighbor angles. Returns: (dominant_symmetry_order, confidence, angle_histogram) """ center = np.mean(unique_2d, axis=0) # Use only central points (avoid edge effects) dists_to_center = np.linalg.norm(unique_2d - center, axis=1) central = unique_2d[dists_to_center < radius] if len(central) < 3: return None, 0.0, None all_nn_angles = [] for pt in central: dists = np.linalg.norm(unique_2d - pt, axis=1) # Exclude self mask = dists > 0.01 neighbor_dists = dists[mask] neighbor_pts = unique_2d[mask] if len(neighbor_dists) < 2: continue # Find k nearest neighbors k = min(k_neighbors, len(neighbor_dists)) idx = np.argsort(neighbor_dists)[:k] nn = neighbor_pts[idx] # Compute angles of neighbor vectors from this point vectors = nn - pt angles = np.arctan2(vectors[:, 1], vectors[:, 0]) angles = np.sort(angles % (2 * np.pi)) # Compute angular differences between consecutive neighbors diffs = np.diff(angles) if len(angles) > 1: # Include wrap-around diffs = np.append(diffs, (2 * np.pi - angles[-1] + angles[0])) all_nn_angles.extend(diffs) if len(all_nn_angles) < 5: return None, 0.0, None all_nn_angles = np.array(all_nn_angles) # Test each symmetry order: what fraction of angular gaps are consistent? scores = {} for n_fold in [2, 3, 4, 6]: expected_gap = 2 * np.pi / n_fold # e.g., 60° for 6-fold # For each observed gap, find the closest multiple of expected_gap # and compute the residual residuals = [] for gap in all_nn_angles: # Which multiple of expected_gap is closest? k_nearest = round(gap / expected_gap) if k_nearest == 0: k_nearest = 1 residual = abs(gap - k_nearest * expected_gap) / expected_gap residuals.append(residual) mean_residual = np.mean(residuals) # Score: 1 - mean_residual (higher = better match) scores[n_fold] = max(0, 1 - mean_residual) # Best match best_n = max(scores, key=scores.get) confidence = scores[best_n] # Map symmetry order to N-wave value # 6-fold rotational symmetry → N=3 (three plane waves at 120° produce 6-fold pattern) # 4-fold rotational symmetry → N=4 # 3-fold rotational symmetry → N=3 # 2-fold → N=2 n_wave_map = {6: 3, 4: 4, 3: 3, 2: 2} return best_n, confidence, scores def compute_nn_symmetry_direct(unique_2d, radius=2.5): """Directly measure rotational symmetry by counting nearest-neighbor angles. More robust method: bin NN angles and look for periodicity. """ center = np.mean(unique_2d, axis=0) dists_to_center = np.linalg.norm(unique_2d - center, axis=1) central = unique_2d[dists_to_center < radius] if len(central) < 5: return None, 0.0, {} # Collect NN distances to find the NN shell all_nn_dists = [] for pt in central: dists = np.linalg.norm(unique_2d - pt, axis=1) dists = dists[dists > 0.01] if len(dists) > 0: all_nn_dists.append(np.min(dists)) if len(all_nn_dists) < 3: return None, 0.0, {} nn_dist = np.median(all_nn_dists) # Adaptive NN shell tolerance: find the gap between 1st and 2nd shells # Collect all pairwise distances to find the shell structure all_pair_dists = [] for pt in central[:min(20, len(central))]: dists = np.linalg.norm(unique_2d - pt, axis=1) dists = np.sort(dists[dists > 0.01])[:8] all_pair_dists.extend(dists) all_pair_dists = np.sort(all_pair_dists) # Find the first significant gap in the sorted distances nn_tol_frac = 0.15 # default 15% tolerance if len(all_pair_dists) > 10: # Look at distances up to 2x the NN distance shell_dists = all_pair_dists[all_pair_dists < nn_dist * 2.5] if len(shell_dists) > 5: gaps = np.diff(shell_dists) # Find the first gap that's >10% of NN distance (shell boundary) big_gaps = np.where(gaps > nn_dist * 0.1)[0] if len(big_gaps) > 0: boundary = shell_dists[big_gaps[0]] + gaps[big_gaps[0]] * 0.5 nn_tol_frac = (boundary / nn_dist) - 1.0 nn_tol_frac = max(0.08, min(nn_tol_frac, 0.25)) nn_tol = nn_dist * nn_tol_frac # Count coordination and collect angles cn_list = [] all_angles = [] for pt in central: dists = np.linalg.norm(unique_2d - pt, axis=1) mask = (dists > 0.01) & (dists < nn_dist + nn_tol) nn_pts = unique_2d[mask] cn_list.append(len(nn_pts)) if len(nn_pts) >= 2: vectors = nn_pts - pt angles = np.arctan2(vectors[:, 1], vectors[:, 0]) angles = np.sort(angles % (2 * np.pi)) # Angular gaps between consecutive neighbors gaps = np.diff(angles) gaps = np.append(gaps, 2 * np.pi - angles[-1] + angles[0]) all_angles.extend(gaps) if len(all_angles) < 3: return None, 0.0, {} all_angles = np.array(all_angles) median_cn = int(np.median(cn_list)) # The median angular gap tells us the symmetry median_gap = np.median(all_angles) median_gap_deg = np.degrees(median_gap) # Score each symmetry order by how well the median gap matches scores = {} for n_fold in [2, 3, 4, 6]: expected_gap_deg = 360.0 / n_fold # Allow matching to the CN-appropriate multiple best_match = float('inf') for mult in [1, 2, 3]: diff = abs(median_gap_deg - expected_gap_deg * mult / max(1, mult)) # Actually just compare directly diff2 = abs(median_gap_deg - 360.0 / n_fold) best_match = min(best_match, diff2) scores[n_fold] = max(0, 1 - best_match / 45.0) # normalize best_symmetry = max(scores, key=scores.get) return best_symmetry, median_cn, scores # ============================================================================ # MAIN ANALYSIS # ============================================================================ def analyze_structure(structure): """Complete analysis for one crystal structure type.""" print(f"\n {'='*60}") print(f" {structure}") print(f" {'='*60}") points_3d, cell = generate_supercell(structure, n_cells=4) directions = get_projection_directions(structure) results = {} for dir_name, normal in directions.items(): # Extract atoms from a SINGLE crystallographic plane # (not a projection of all atoms — that conflates layers) plane_2d, d_spacing, n_layers = extract_single_plane( points_3d, normal / np.linalg.norm(normal), tol=0.05 ) if plane_2d is None or len(plane_2d) < 3: print(f" {dir_name:10s}: insufficient atoms on single plane") results[dir_name] = { 'density': 0, 'symmetry': None, 'median_cn': 0, 'scores': {}, 'n_unique': 0, 'unique_2d': np.zeros((0, 2)), 'd_spacing': d_spacing, 'n_layers': n_layers, } continue # Density of this single plane density = compute_planar_density(plane_2d, radius=2.5) # Symmetry of the atomic arrangement on this plane sym_order, median_cn, scores = compute_nn_symmetry_direct( plane_2d, radius=2.0 ) results[dir_name] = { 'density': density, 'symmetry': sym_order, 'median_cn': median_cn, 'scores': scores, 'n_unique': len(plane_2d), 'unique_2d': plane_2d, 'd_spacing': d_spacing, 'n_layers': n_layers, } sym_str = f"{sym_order}-fold" if sym_order else "N/A" scores_str = ", ".join(f"N={k}:{v:.2f}" for k, v in sorted(scores.items())) print(f" {dir_name:10s}: density={density:.3f} CN_2D={median_cn} " f"symmetry={sym_str} n_atoms={len(plane_2d)} " f"d_spacing={d_spacing:.3f} [{scores_str}]") # Find densest plane densest_dir = max(results, key=lambda d: results[d]['density']) densest_sym = results[densest_dir]['symmetry'] # Map rotational symmetry to N-wave value # 6-fold → N=3 (three waves at 120° create 6-fold pattern) # 4-fold → N=4 # 3-fold → N=3 # 2-fold → N=2 sym_to_N = {6: 3, 4: 4, 3: 3, 2: 2} blind_N = sym_to_N.get(densest_sym, None) model_N = MODEL_MAPPING.get(structure) match = "MATCH" if blind_N == model_N else "MISMATCH" print(f"\n Densest plane: {densest_dir} (density={results[densest_dir]['density']:.3f})") print(f" Densest plane symmetry: {densest_sym}-fold → blind N = {blind_N}") print(f" Model assignment: N = {model_N}") print(f" *** {match} ***") return { 'structure': structure, 'projections': results, 'densest_dir': densest_dir, 'densest_sym': densest_sym, 'blind_N': blind_N, 'model_N': model_N, 'match': blind_N == model_N, } def generate_projection_plots(all_results): """Generate visualization of 2D projections for each structure.""" n_structures = len(all_results) fig, axes = plt.subplots(n_structures, 3, figsize=(15, 4 * n_structures)) fig.suptitle("TLT-006: Blind N-Fold Assignment — 2D Crystal Projections\n" "Green border = densest plane (blind N source) | " "Title shows detected symmetry", fontsize=13, fontweight='bold') for row, res in enumerate(all_results): directions = list(res['projections'].keys()) for col, dir_name in enumerate(directions[:3]): ax = axes[row, col] if n_structures > 1 else axes[col] proj = res['projections'][dir_name] pts = proj['unique_2d'] # Plot points ax.scatter(pts[:, 0], pts[:, 1], s=15, c='navy', alpha=0.7, edgecolors='black', linewidths=0.3) # Highlight if densest is_densest = (dir_name == res['densest_dir']) if is_densest: for spine in ax.spines.values(): spine.set_color('green') spine.set_linewidth(3) sym = proj['symmetry'] cn = proj['median_cn'] dens = proj['density'] sym_str = f"{sym}-fold" if sym else "N/A" match_str = "" if is_densest: match_str = " [DENSEST]" if res['match']: match_str += " ✓" else: match_str += " ✗" ax.set_title(f"{res['structure']} {dir_name}\n" f"CN={cn} {sym_str} d={dens:.2f}{match_str}", fontsize=9, fontweight='bold' if is_densest else 'normal', color='green' if is_densest and res['match'] else 'red' if is_densest and not res['match'] else 'black') ax.set_aspect('equal') ax.set_xlim(-3, 3) ax.set_ylim(-3, 3) ax.grid(True, alpha=0.2) ax.tick_params(labelsize=6) plt.tight_layout(rect=[0, 0, 1, 0.95]) plt.savefig(OUTPUT_DIR / "blind_N_projections.png", dpi=150, bbox_inches='tight') plt.close() def generate_summary_table(all_results): """Generate summary comparison table.""" fig, ax = plt.subplots(figsize=(12, 5)) ax.axis('off') headers = ['Structure', 'Densest\nPlane', 'Density', 'CN_2D', 'Symmetry', 'Blind N', 'Model N', 'Match?'] cell_data = [] for r in all_results: densest = r['projections'][r['densest_dir']] cell_data.append([ r['structure'], r['densest_dir'], f"{densest['density']:.3f}", str(densest['median_cn']), f"{r['densest_sym']}-fold", str(r['blind_N']), str(r['model_N']), 'YES' if r['match'] else 'NO', ]) colors = [['lightgreen' if row[-1] == 'YES' else '#ffcccc'] * len(headers) for row in cell_data] table = ax.table(cellText=cell_data, colLabels=headers, cellColours=colors, loc='center', cellLoc='center') table.auto_set_font_size(False) table.set_fontsize(10) table.scale(1, 1.8) ax.set_title("TLT-006: Blind N-Fold Assignment — Summary\n" "Does the densest crystallographic plane's symmetry match " "the model's N assignment?", fontsize=12, fontweight='bold', pad=20) plt.savefig(OUTPUT_DIR / "blind_N_summary.png", dpi=150, bbox_inches='tight') plt.close() def run_analysis(): print("=" * 72) print("TLT-006: BLIND N-FOLD ASSIGNMENT VALIDATION") print("=" * 72) print(""" PURPOSE: Validate the N-fold symmetry mapping used in TLT-004 without relying on the mapping itself. Instead, examine each crystal structure's 2D projections along all high-symmetry directions and let the DATA determine which N best matches. METHOD: 1. Generate 3D supercell for each standard crystal type 2. Project along [100], [110], [111] (cubic) or [0001], [10-10], [11-20] (hex) 3. Compute 2D planar density and nearest-neighbor angular symmetry 4. The densest plane's symmetry determines the "blind N" 5. Compare blind N against model N CRITERION: Densest crystallographic plane (objective, physically motivated) - Densest plane determines cleavage, growth, slip, and surface stability - Not a free parameter — uniquely determined by the crystal structure MODEL ASSIGNMENTS BEING TESTED: FCC → N=3, HCP → N=3, BCC → N=4, diamond → N=3, SC → N=4, BCT → N=4, rhombo → N=3 """) structures = ['FCC', 'BCC', 'HCP', 'diamond', 'SC', 'BCT', 'rhombo'] all_results = [] for struct in structures: result = analyze_structure(struct) all_results.append(result) # Summary n_match = sum(1 for r in all_results if r['match']) n_total = len(all_results) print("\n" + "=" * 72) print("SUMMARY") print("=" * 72) print(f"\n Concordance: {n_match}/{n_total} ({100*n_match/n_total:.0f}%)") print(f"\n {'Structure':10s} {'Densest':10s} {'Blind N':8s} {'Model N':8s} {'Match':6s}") print(f" {'-'*42}") for r in all_results: match_str = "YES" if r['match'] else "*** NO ***" print(f" {r['structure']:10s} {r['densest_dir']:10s} " f"N={r['blind_N']!s:5s} N={r['model_N']!s:5s} {match_str}") # Detail on mismatches mismatches = [r for r in all_results if not r['match']] if mismatches: print(f"\n MISMATCHES ({len(mismatches)}):") for r in mismatches: print(f"\n {r['structure']}:") print(f" Model says N={r['model_N']}") print(f" Densest plane {r['densest_dir']} has " f"{r['densest_sym']}-fold symmetry → N={r['blind_N']}") # Show alternative planes for dir_name, proj in r['projections'].items(): sym_to_N = {6: 3, 4: 4, 3: 3, 2: 2} proj_N = sym_to_N.get(proj['symmetry']) match_model = " ← matches model" if proj_N == r['model_N'] else "" print(f" {dir_name}: {proj['symmetry']}-fold (N={proj_N}) " f"density={proj['density']:.3f}{match_model}") else: print("\n *** ALL ASSIGNMENTS MATCH ***") print(" The model's N-fold mapping is empirically validated by the") print(" densest-plane criterion for all tested crystal structures.") # Interpretation print(f""" INTERPRETATION: The blind test assigns N based on the densest crystallographic plane's rotational symmetry. This is physically motivated: the densest plane is the most stable and determines the dominant geometry. Concordance rate: {n_match}/{n_total} ({100*n_match/n_total:.0f}%) {"The model's mapping is SUPPORTED by the densest-plane criterion." if n_match == n_total else "Some assignments DISAGREE with the densest-plane criterion."} {"" if n_match == n_total else "The model mapping for mismatched structures should be flagged as choices rather than derivations."} IMPORTANT CAVEAT: Even for matched structures, the mapping choice (densest plane) is one of several physically reasonable criteria. An alternative criterion (e.g., highest-symmetry axis, most common plane family) might yield different results. The test shows that the model's mapping is CONSISTENT WITH a specific physical rule, not that it is the UNIQUE correct mapping. AUDIT STATUS: This test provides empirical support for the mapping scored 2/10 by Gemini and 6/10 by Grok. The densest-plane rule is objective and reproducible. If {n_match}/{n_total} structures match, the mapping is not arbitrary — it corresponds to a physical property, even if not derived from first principles. """) # Generate plots generate_projection_plots(all_results) generate_summary_table(all_results) print(f" Plots saved to {OUTPUT_DIR}/") # Write data file write_data_file(all_results) return all_results def write_data_file(all_results): """Write detailed data file.""" filepath = OUTPUT_DIR / "TLT-006_blind_N_data.txt" with open(filepath, 'w') as f: f.write("=" * 72 + "\n") f.write("TLT-006: BLIND N-FOLD ASSIGNMENT VALIDATION — DATA\n") f.write("=" * 72 + "\n") f.write(f"Date: 2026-03-13\n\n") f.write("METHOD: Densest crystallographic plane symmetry\n") f.write("CRITERION: Objective — densest plane is uniquely determined\n\n") for r in all_results: f.write(f"{r['structure']}:\n") f.write(f" Model N = {r['model_N']}\n") f.write(f" Blind N = {r['blind_N']} (from {r['densest_dir']})\n") f.write(f" Match: {'YES' if r['match'] else 'NO'}\n") f.write(f" Projections:\n") for dir_name, proj in r['projections'].items(): densest = " [DENSEST]" if dir_name == r['densest_dir'] else "" f.write(f" {dir_name}: density={proj['density']:.4f} " f"CN_2D={proj['median_cn']} " f"symmetry={proj['symmetry']}-fold " f"(scores: {proj['scores']}){densest}\n") f.write("\n") n_match = sum(1 for r in all_results if r['match']) f.write(f"\nConcordance: {n_match}/{len(all_results)}\n") print(f" Data file written: {filepath}") if __name__ == "__main__": run_analysis() print("\n" + "=" * 72) print("TLT-006 COMPLETE") print("=" * 72)