逐步聚合物构建¶

import molpy as mp
import numpy as np
from pathlib import Path
from molpy.reacter import find_neighbors, find_port
from molpy.typifier import OplsTypifier

output_dir = Path("02_output")
output_dir.mkdir(exist_ok=True)

ff = mp.io.read_xml_forcefield("oplsaa.xml")
typifier = OplsTypifier(ff, strict_typing=False)

MONOMER_BIGSMILES = "{[][<]OCCOCCOCCO[>][]}"


def make_eo_monomer():
    m = mp.parser.parse_monomer(MONOMER_BIGSMILES)
    m = mp.adapter.RDKitAdapter(m).generate_3d(add_hydrogens=True, optimize=True)
    m = m.get_topo(gen_angle=True, gen_dihe=True)
    m = typifier.typify(m)
    return m


mon_a = make_eo_monomer()
mon_b = make_eo_monomer()
mon_b.move([10.0, 0.0, 0.0])  # 平移到远离 merged 后坐标重叠的位置

print(f"单体: {len(mon_a.atoms)} 个原子, {len(mon_a.bonds)} 根键")
print(f"      {len(mon_a.angles)} 个键角, {len(mon_a.dihedrals)} 个二面角")
print(f"未分类: {sum(1 for a in mon_a.atoms if a.get('type') is None)} 个原子")

def save_step(name: str, chain, ff, output_dir: Path) -> None:
    frame = chain.to_frame()
    atoms = frame["atoms"]
    atoms["mol_id"] = np.ones(atoms.nrows, dtype=int)
    coords = np.column_stack([atoms["x"], atoms["y"], atoms["z"]])
    lo = coords.min(axis=0) - 5.0
    hi = coords.max(axis=0) + 5.0
    frame.box = mp.Box(matrix=hi - lo, origin=lo)
    mp.io.write_lammps_system(output_dir / name, frame, ff)
    print(f"  已保存: {output_dir / name / (name + '.data')}")


save_step("peo1", mon_a, ff, output_dir)

# mon_a 的右端口和 mon_b 的左端口
port_r = find_port(mon_a, ">")
port_l = find_port(mon_b, "<")

# 锚定碳是右端口氧的碳邻居
anchor_C = [a for a in find_neighbors(mon_a, port_r) if a.get("element") == "C"][0]

# 离去基团 –OH：端口氧及其单个氢
leaving_OH_oxygen = port_r
leaving_OH_hydrogen = find_neighbors(mon_a, leaving_OH_oxygen, element="H")[0]

# 离去 H：mon_b 左侧端口上的一个氢
leaving_H_b = find_neighbors(mon_b, port_l, element="H")[0]

print(f"新键将形成:  {anchor_C.get('element')} — {port_l.get('element')}")
print(
    f"从 mon_a 离去:  {leaving_OH_oxygen.get('element')} + {leaving_OH_hydrogen.get('element')}  (= OH)"
)
print(f"从 mon_b 离去:  {leaving_H_b.get('element')}  (= H)")
print(f"净脱除:         H₂O")

# 1. 将两个单体合并为单个 Atomistic（两个单体的原子和键都合到一起）
dimer = mon_a.merge(mon_b)
print(
    f"合并后: {len(dimer.atoms)} 个原子  (= {len(mon_a.atoms)} + {len(mon_b.atoms)})"
)

# 2. 在锚定碳和 mon_b 的端口氧之间形成新 C–O 键
new_bond = dimer.def_bond(anchor_C, port_l)
print(
    f"新键: {new_bond.endpoints[0].get('element')}—{new_bond.endpoints[1].get('element')}"
)

# 3. 移除三个离去原子；它们关联的键会自动删除
dimer.del_atom(leaving_OH_oxygen, leaving_OH_hydrogen, leaving_H_b)
print(
    f"移除后: {len(dimer.atoms)} 个原子  (= {len(dimer.atoms)} = {len(mon_a.atoms) + len(mon_b.atoms)} − 3)"
)

# 检查修复前连接区域中未解析的内容
print(
    f"get_topo 前未分类的键:      {sum(1 for b in dimer.bonds if b.get('type') is None)}"
)
print(
    f"get_topo 前未分类的键角:     {sum(1 for a in dimer.angles if a.get('type') is None)}"
)
print(
    f"get_topo 前未分类的二面角:  {sum(1 for d in dimer.dihedrals if d.get('type') is None)}"
)

angles_before = len(dimer.angles)
dihe_before = len(dimer.dihedrals)
dimer = dimer.get_topo(gen_angle=True, gen_dihe=True)

print(
    f"\n键角:    {angles_before} → {len(dimer.angles)}  (+{len(dimer.angles) - angles_before} 跨连接)"
)
print(
    f"二面角: {dihe_before} → {len(dimer.dihedrals)}  (+{len(dimer.dihedrals) - dihe_before} 跨连接)"
)

dimer = typifier.typify(dimer)

print(f"\n未分类的键:      {sum(1 for b in dimer.bonds if b.get('type') is None)}")
print(f"未分类的键角:     {sum(1 for a in dimer.angles if a.get('type') is None)}")
print(f"未分类的二面角:  {sum(1 for d in dimer.dihedrals if d.get('type') is None)}")

save_step("peo2", dimer, ff, output_dir)

chain = make_eo_monomer()
save_step("peo1", chain, ff, output_dir)

for i in range(1, 5):
    unit = make_eo_monomer()
    unit.move([10.0 * i, 0.0, 0.0])

    # 识别当前链端和传入单元上的反应原子
    p_r = find_port(chain, ">")
    p_l = find_port(unit, "<")
    anc = [a for a in find_neighbors(chain, p_r) if a.get("element") == "C"][0]
    l_O = p_r
    l_H1 = find_neighbors(chain, l_O, element="H")[0]
    l_H2 = find_neighbors(unit, p_l, element="H")[0]

    # 偶合
    chain = chain.merge(unit)
    chain.def_bond(anc, p_l)
    chain.del_atom(l_O, l_H1, l_H2)
    chain = chain.get_topo(gen_angle=True, gen_dihe=True)
    chain = typifier.typify(chain)

    name = f"peo{i + 1}"
    save_step(name, chain, ff, output_dir)
    print(f"{name}: {len(chain.atoms)} 个原子, {len(chain.bonds)} 根键")

from molpy.core.atomistic import Atom, Atomistic
from molpy.reacter import (
    Reacter,
    form_single_bond,
    select_hydrogens,
    select_neighbor,
    select_self,
)


def select_hydroxyl_group(struct: Atomistic, reaction_site: Atom) -> list[Atom]:
    """返回 [O, H] — 锚定碳上的羟基离去基团。"""
    for neighbor in find_neighbors(struct, reaction_site):
        if neighbor.get("element") != "O":
            continue
        h_neighbors = find_neighbors(struct, neighbor, element="H")
        if h_neighbors:
            return [neighbor, h_neighbors[0]]
    raise ValueError("反应位点附近未找到羟基")


rxn = Reacter(
    name="dehydration",
    anchor_selector_left=select_neighbor("C"),  # 等同于：找到端口 > 旁边的 C
    anchor_selector_right=select_self,  # 等同于：端口 < 的 O 本身
    leaving_selector_left=select_hydroxyl_group,  # 等同于：l_O + l_H1
    leaving_selector_right=select_hydrogens(1),  # 等同于：l_H2
    bond_former=form_single_bond,  # 等同于：def_bond(anc, p_l)
)

chain = make_eo_monomer()
save_step("peo1_rxn", chain, ff, output_dir)

for i in range(1, 5):
    unit = make_eo_monomer()
    unit.move([10.0 * i, 0.0, 0.0])

    result = rxn.run(
        left=chain,
        right=unit,
        port_atom_L=find_port(chain, ">"),
        port_atom_R=find_port(unit, "<"),
        compute_topology=True,
    )
    chain = result.product
    chain = chain.get_topo(gen_angle=True, gen_dihe=True)
    chain = typifier.typify(chain)

    name = f"peo{i + 1}_rxn"
    save_step(name, chain, ff, output_dir)
    print(f"{name}: {len(chain.atoms)} 个原子")

for a in mon_a.atoms:
    if a.get("port"):
        print(f"  element={a.get('element')}  port={a.get('port')}")

p_r = find_port(mon_a, ">")
anc = [a for a in find_neighbors(mon_a, p_r) if a.get("element") == "C"][0]
for nb in find_neighbors(mon_a, anc):
    print(f"  element={nb.get('element')}")

chain = chain.merge(unit)
chain.def_bond(anc, p_l)
chain.del_atom(l_O, l_H1, l_H2)
chain.get_topo(gen_angle=True, gen_dihe=True)   # 第一步
chain = typifier.typify(chain)                  # 第二步

frame["atoms"]["mol_id"] = np.ones(frame["atoms"].nrows, dtype=int)