拓扑驱动的 CGSmiles 组装¶

from molpy.parser import parse_cgsmiles

expressions = {
    "linear": "{[#EO2]|4[#PS]}",
    "ring": "{[#EO2]1[#PS][#EO2][#PS][#EO2]1}",
    "branch": "{[#PS][#EO3]([#PS])[#PS]}",
}

for name, expr in expressions.items():
    ir = parse_cgsmiles(expr)
    labels = [node.label for node in ir.base_graph.nodes]
    print(f"{name}: nodes={len(ir.base_graph.nodes)}, labels={labels}")

import molpy as mp
from molpy.typifier import OplsTypifier

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

BIGSMILES = {
    "EO2": "{[][$]OCCO[$][]}",
    "EO3": "{[][$]OCC(CO[$])(CO[$])[]}",
    "PS": "{[][$]OCC(c1ccccc1)CO[$][]}",
}


def build_monomer(bigsmiles, typifier):
    monomer = mp.parser.parse_monomer(bigsmiles)
    monomer = mp.adapter.RDKitAdapter(monomer).generate_3d(add_hydrogens=True, optimize=False)
    monomer = monomer.get_topo(gen_angle=True, gen_dihe=True)
    monomer = typifier.typify(monomer)
    return monomer


library = {label: build_monomer(bs, typifier) for label, bs in BIGSMILES.items()}

for label, mon in library.items():
    ports = [a.get("port") for a in mon.atoms if a.get("port")]
    print(f"{label}: atoms={len(mon.atoms)}, ports={ports}")

from molpy.core.atomistic import Atom, Atomistic
from molpy.builder.polymer import (
    Connector,
    CovalentSeparator,
    LinearOrienter,
    Placer,
    PolymerBuilder,
)
from molpy.reacter import (
    Reacter,
    find_neighbors,
    form_single_bond,
    select_neighbor,
    select_self,
)


def select_hydroxyl_group(struct: Atomistic, site: Atom) -> list[Atom]:
    for nb in find_neighbors(struct, site):
        if nb.get("element") != "O":
            continue
        hs = [a for a in find_neighbors(struct, nb, element="H")]
        if hs:
            return [nb, hs[0]]
    raise ValueError("No hydroxyl group found")


def select_one_hydrogen(struct: Atomistic, site: Atom) -> list[Atom]:
    hs = [a for a in find_neighbors(struct, site, element="H")]
    if not hs:
        raise ValueError("No hydrogen found")
    return [hs[0]]


rxn = Reacter(
    name="dehydration",
    anchor_selector_left=select_neighbor("C"),
    anchor_selector_right=select_self,
    leaving_selector_left=select_hydroxyl_group,
    leaving_selector_right=select_one_hydrogen,
    bond_former=form_single_bond,
)

rules = {(l, r): ("$", "$") for l in library for r in library}
connector = Connector(port_map=rules, reacter=rxn)
placer = Placer(
    separator=CovalentSeparator(buffer=-0.1),
    orienter=LinearOrienter(),
)

builder = PolymerBuilder(
    library=library,
    connector=connector,
    placer=placer,
    typifier=typifier,
)

for name, expr in expressions.items():
    result = builder.build(expr)
    polymer = result.polymer
    print(f"{name}: atoms={len(polymer.atoms)}, steps={result.total_steps}")

import numpy as np
from pathlib import Path

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

for name, expr in expressions.items():
    result = builder.build(expr)
    typed_polymer = typifier.typify(result.polymer)
    frame = typed_polymer.to_frame()
    atoms = frame["atoms"]
    if "mol_id" not in 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)