import numpy as np
import matplotlib.pyplot as plt
import itertools
import random
from matplotlib.animation import FuncAnimation
# ====== 常量 ======
DIRECTIONS = [
np.array([1.0, 0.0]),
np.array([0.5, np.sqrt(3) / 2]),
np.array([-0.5, np.sqrt(3) / 2])
]
ROD_LEN = 0.8
BOND_GAP = 0.08
ATOM_CLR = "#2e6bd9" # 自由原子
LOCKED_CLR = "green" # 已成键原子
ENDPOINT_CLR = "orange" # 链条端点
BOND_COLOR = {"A": "red", "B": "yellow"} # 区分 A / B
# ====== 三角格点 ======
def tri_lattice_points(L):
pts = []
for j in range(L):
for i in range(L):
x = i + 0.5 * (j % 2)
y = j * np.sqrt(3) / 2
pts.append((x, y))
return np.asarray(pts)
def lattice_neighbors(L, index):
pts = tri_lattice_points(L)
c = pts[index]
neighbors = []
for j, p in enumerate(pts):
if np.isclose(np.linalg.norm(c - p), 1.0, atol=1e-6):
neighbors.append(j)
return neighbors
# ====== 判定函数 ======
def equilateral_side1(c1, c2, c3, tol=1e-3):
d12 = np.linalg.norm(c1 - c2)
d23 = np.linalg.norm(c2 - c3)
d31 = np.linalg.norm(c3 - c1)
return (np.isclose(d12, 1.0, atol=tol) and
np.isclose(d23, 1.0, atol=tol) and
np.isclose(d31, 1.0, atol=tol))
def is_collinear(c1, c2, c3):
if np.isclose(c1[1], c2[1]) and np.isclose(c2[1], c3[1]):
return True
vec1 = c2 - c1
vec2 = c3 - c1
cross = vec1[0] * vec2[1] - vec1[1] * vec2[0]
return np.isclose(cross, 0, atol=1e-6)
def triangle_type(c1, d1, c2, d2, c3, d3):
if is_collinear(c1, c2, c3):
return None
dirs = {d1, d2, d3}
if len(dirs) == 1:
return "A"
if len(dirs) == 3:
edges = [c2 - c1, c3 - c2, c1 - c3]
edges = [e / np.linalg.norm(e) for e in edges]
centers = [(c1, d1), (c2, d2), (c3, d3)]
for center, d in centers:
v = DIRECTIONS[d] / np.linalg.norm(DIRECTIONS[d])
if not any(abs(np.dot(v, e)) > 0.99 for e in edges):
return None
for other, _ in centers:
if np.allclose(other, center):
continue
vec = other - center
cross = abs(np.cross(v, vec / np.linalg.norm(vec)))
if cross < 1e-6:
return None
return "B"
return None
# ====== 原子绘制 ======
def atom_endpoints(center, dir_idx, rod_len=ROD_LEN, inner_gap=BOND_GAP):
v = DIRECTIONS[dir_idx] / np.linalg.norm(DIRECTIONS[dir_idx])
half = rod_len / 2.0 - inner_gap
return center - half * v, center + half * v
def closest_side_connection(centerA, dirA, centerB, dirB):
a1, a2 = atom_endpoints(centerA, dirA)
b1, b2 = atom_endpoints(centerB, dirB)
candidates = [
(np.linalg.norm(a1 - b1), a1, b1),
(np.linalg.norm(a1 - b2), a1, b2),
(np.linalg.norm(a2 - b1), a2, b1),
(np.linalg.norm(a2 - b2), a2, b2),
]
_, pa, pb = min(candidates, key=lambda x: x[0])
return pa, pb
# ====== 仿真类 ======
class Simulation:
def __init__(self, L=6, N=12, seed=42):
np.random.seed(seed)
self.L = L
self.points = tri_lattice_points(L)
idxs = np.random.choice(len(self.points), size=N, replace=False)
self.atoms = [(idx, np.random.randint(0, 3)) for idx in idxs]
self.locked = [False] * N
self.used_in_B = [False] * N
self.triangles = []
self.has_B = False
self.locked_edges = set()
self.chain_atoms = set()
self.chain_triangles = []
self.chain_endpoints = set()
def step(self):
new_tris = []
coords = [self.points[idx] for idx, _ in self.atoms]
dirs = [d for _, d in self.atoms]
for (i, j, k) in itertools.combinations(range(len(self.atoms)), 3):
# 只有 B 类原子才会被 used_in_B 标记，因此它们不会参与后续成键
if self.used_in_B[i] or self.used_in_B[j] or self.used_in_B[k]:
continue
edges = [(min(i, j), max(i, j)), (min(j, k), max(j, k)), (min(k, i), max(k, i))]
if any(edge in self.locked_edges for edge in edges):
continue
c1, c2, c3 = coords[i], coords[j], coords[k]
if equilateral_side1(c1, c2, c3):
t = triangle_type(c1, dirs[i], c2, dirs[j], c3, dirs[k])
if t == "B" and not self.has_B:
new_tris.append((i, j, k, "B"))
self.has_B = True
for idx in (i, j, k):
self.used_in_B[idx] = True
self.locked[idx] = True
self.locked_edges.update(edges)
elif t == "A" and self.has_B:
tri_set = {i, j, k}
if not self.chain_triangles:
# 第一条A链，三原子都没被A链用过即可
if all(idx not in self.chain_atoms for idx in tri_set):
new_tris.append((i, j, k, "A"))
self.chain_triangles.append((i, j, k))
self.chain_atoms.update(tri_set)
self.chain_endpoints = set(tri_set)
for idx in (i, j, k):
self.locked[idx] = True
self.locked_edges.update(edges)
else:
# 只允许新A型三角形与A链端点和一个新原子组成
endpoints = self.chain_endpoints
inter = tri_set & endpoints
if len(inter) == 2 and len(tri_set - self.chain_atoms) == 1:
new_tris.append((i, j, k, "A"))
self.chain_triangles.append((i, j, k))
self.chain_atoms.update(tri_set)
new_atom = list(tri_set - self.chain_atoms)[0]
covered = list(inter)
# 计算链条中心
chain_centers = [self.points[self.atoms[idx][0]] for idx in self.chain_atoms]
chain_center = np.mean(chain_centers, axis=0)
# 计算两个被覆盖端点到链中心的距离
dists = [np.linalg.norm(self.points[self.atoms[c][0]] - chain_center) for c in covered]
# 靠近链中心的端点变为非端点，远离中心的端点保留
min_idx = np.argmin(dists)
# 移除靠近中心的端点
self.chain_endpoints.discard(covered[min_idx])
# 新增新原子为端点
self.chain_endpoints.add(new_atom)
for idx in (i, j, k):
self.locked[idx] = True
self.locked_edges.update(edges)
# ====== 强制生长A链 ======
# 如果A链已存在且有端点，强制将一个自由原子移动到端点附近并尝试生成A型三角形
if self.chain_endpoints and len(self.chain_atoms) < len(self.atoms):
# 找到所有未锁定的自由原子
free_idxs = [idx for idx, locked in enumerate(self.locked) if not locked and idx not in self.chain_atoms]
if free_idxs and self.chain_endpoints:
# 任选一个端点和一个自由原子
endpoint = next(iter(self.chain_endpoints))
free_idx = free_idxs[0]
# 让自由原子移动到端点最近的邻居格点，并方向与端点原子一致
endpoint_site, endpoint_dir = self.atoms[endpoint]
endpoint_pos = self.points[endpoint_site]
neigh = lattice_neighbors(self.L, endpoint_site)
# 只选未被占用的邻居
occupied_sites = {self.atoms[idx][0] for idx, locked in enumerate(self.locked) if locked}
occupied_sites |= {self.atoms[idx][0] for idx in self.chain_atoms}
free_neigh = [n for n in neigh if n not in occupied_sites]
if free_neigh:
new_site = free_neigh[0]
self.atoms[free_idx] = (new_site, endpoint_dir)
# 立即尝试生成新的A型三角形（下次step会自动检测）
self.triangles.extend(new_tris)
occupied_sites = {self.atoms[idx][0] for idx, locked in enumerate(self.locked) if locked}
for idx, locked in enumerate(self.locked):
if not locked:
site, d = self.atoms[idx]
neigh = lattice_neighbors(self.L, site)
free_neigh = [n for n in neigh if n not in occupied_sites]
if free_neigh:
site = random.choice(free_neigh)
d = random.randint(0, 2)
self.atoms[idx] = (site, d)
occupied_sites.add(site)
return all(self.locked)
# ====== 绘制函数 ======
def plot_frame(sim, ax):
ax.clear()
coords = [sim.points[idx] for idx, _ in sim.atoms]
dirs = [d for _, d in sim.atoms]
atom_counts = {a: 0 for a in sim.chain_atoms}
for tri in sim.chain_triangles:
for v in tri:
atom_counts[v] += 1
endpoints = {a for a, c in atom_counts.items() if c == 1}
for (center, dir_idx), locked, idx in zip(zip(coords, dirs), sim.locked, range(len(sim.atoms))):
p1, p2 = atom_endpoints(center, dir_idx)
if idx in endpoints:
color = ENDPOINT_CLR
else:
color = LOCKED_CLR if locked else ATOM_CLR
ax.plot([p1[0], p2[0]], [p1[1], p2[1]],
color=color, lw=6, solid_capstyle="round")
for (i, j, k, t) in sim.triangles:
edges = [(i, j), (j, k), (k, i)]
for a, b in edges:
pa, pb = closest_side_connection(coords[a], dirs[a], coords[b], dirs[b])
ax.plot([pa[0], pb[0]], [pa[1], pb[1]], color=BOND_COLOR[t], lw=2)
# 固定画幅为 9×9
ax.set_xlim(-0.5, 8.5)
ax.set_ylim(-0.5, 8.5)
ax.set_aspect("equal")
# ====== 动画 ======
def run_animation():
sim = Simulation(L=6, N=12, seed=1)
fig, ax = plt.subplots(figsize=(6, 6))
def update(frame):
done = sim.step()
plot_frame(sim, ax)
if done:
ani.event_source.stop()
ani = FuncAnimation(fig, update, frames=500, interval=200)
plt.show()
# ====== 运行动画 ======
run_animation()我现在遇到了问题，A型键无法生长
最新发布
本文详细介绍了Python中的几个关键函数，包括`itertools.combinations()`用于生成集合和列表的无重复组合，`combinations_with_replacement()`允许元素自我组合，`zip()`进行数据压缩与解压缩，`format()`进行字符串格式化，`enumerate()`生成索引序列，以及`np.linalg.norm()`计算范数。
扫一扫
618
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
