交联网络¶

from pathlib import Path
import numpy as np
import molpy as mp
from molpy.core.atomistic import Atom, Atomistic
from molpy.core.element import Element
from molpy.typifier import OplsTypifier

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

eo2 = mp.parser.parse_monomer("{[][$]OCCOCCOCCO[$][]}")
eo2 = mp.adapter.RDKitAdapter(eo2).generate_3d(add_hydrogens=True, optimize=True)
eo2 = eo2.get_topo(gen_angle=True, gen_dihe=True)
for idx, atom in enumerate(eo2.atoms, start=1):
    atom["id"] = idx
eo2 = typifier.typify(eo2)

eo3 = mp.parser.parse_monomer("{[]C(COCCO[$])(COCCO[$])COCCO[$][]}")
eo3 = mp.adapter.RDKitAdapter(eo3).generate_3d(add_hydrogens=True, optimize=True)
eo3 = eo3.get_topo(gen_angle=True, gen_dihe=True)
for idx, atom in enumerate(eo3.atoms, start=1):
    atom["id"] = idx
eo3 = typifier.typify(eo3)

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

from molpy.reacter import (
    BondReactReacter,
    find_neighbors,
    find_port,
    form_single_bond,
    select_hydrogens,
    select_hydroxyl_group,
    select_neighbor,
    select_self,
)

reacter = BondReactReacter(
    name="rxn1",
    anchor_selector_left=select_neighbor("C"),
    anchor_selector_right=select_self,
    leaving_selector_left=select_hydroxyl_group,
    leaving_selector_right=select_hydrogens(1),
    bond_former=form_single_bond,
    radius=4,
)

left = eo2.copy()
right = eo2.copy()
result = reacter.run(
    left=left,
    right=right,
    port_atom_L=find_port(left, "$"),
    port_atom_R=find_port(right, "$"),
    compute_topology=True,
    typifier=typifier,
)
template = result.template

print(f"pre:  {len(template.pre.atoms)} atoms")
print(f"post: {len(template.post.atoms)} atoms")

from molpy.pack import InsideBoxConstraint, Packmol

N_EO2, N_EO3 = 27, 9
TARGET_DENSITY = 1.1  # g/cm³ (无定形 PEO ≈ 1.1–1.2)

total_mass_g = (
    N_EO2 * sum(Element(a.get("element")).mass for a in eo2.atoms)
    + N_EO3 * sum(Element(a.get("element")).mass for a in eo3.atoms)
) / 6.022e23
box_length = ((total_mass_g / TARGET_DENSITY) * 1e24) ** (1 / 3)

packer = Packmol(workdir=Path("04_output/packmol"))
constraint = InsideBoxConstraint(
    length=np.array([box_length] * 3),
    origin=np.zeros(3),
)
packer.def_target(eo2.to_frame(), number=N_EO2, constraint=constraint)
packer.def_target(eo3.to_frame(), number=N_EO3, constraint=constraint)

packed = packer(max_steps=20000, seed=42)
packed.box = mp.Box.cubic(length=box_length)
print(f"packed: {packed['atoms'].nrows} atoms in {box_length:.1f} Å box")

from pathlib import Path

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

mp.io.write_lammps_bond_react_system(
    output_dir,
    packed,
    ff,
    templates={"rxn1": template},
)

lammps_script = """\
# ====================================================================
# PEO 交联网络 – OPLS-AA / fix bond/react
# ====================================================================
units           real
atom_style      full
boundary        p p p

read_data       04_output.data
include         04_output.ff

# -- 长程静电 --
kspace_style    pppm 1.0e-4

# -- OPLS-AA 1-4 缩放 --
special_bonds   lj 0.0 0.0 0.5 coul 0.0 0.0 0.5

neighbor        2.0 bin
neigh_modify    every 1 delay 0 check yes

# ====================================================================
# 阶段 1：能量最小化（消除堆积重叠）
# ====================================================================
thermo          100
thermo_style    custom step temp pe ke etotal press vol density
min_style       cg
minimize        1.0e-4 1.0e-6 2000 20000

reset_timestep  0

# ====================================================================
# 注意：运行步数已缩短（每阶段仅几步），以便本文档页面快速渲染；
# 实际交联运行请使用 5 / 5 / 10+ ps。
# 阶段 2：NVT 平衡（300 K）
# ====================================================================
variable step equal "step"
variable T    equal "temp"
variable rho  equal "density"
variable rxn1 equal "f_rxns[1]"

fix out all print 100 "${step} ${T} ${rho}" &
    file thermo_rxn.dat screen no &
    title "# step temp density"
velocity        all create 300.0 12345 dist gaussian

fix             nvt_eq all nvt temp 300.0 300.0 100.0
timestep        1.0
run             200
unfix           nvt_eq

# ====================================================================
# 阶段 3：NPT 平衡（300 K, 1 atm）
# ====================================================================
fix             equil all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
run             200
run             200
run             100
unfix           equil

# ====================================================================
# 阶段 4：使用 fix bond/react 进行反应性 MD（300 K NPT）
# ====================================================================
molecule        rxn1_pre rxn1_pre.mol
molecule        rxn1_post rxn1_post.mol

# bond/react: 每一步尝试，5 Å 截断，仅分子间反应
fix             rxns all bond/react stabilization yes npt_grp 0.03 &
                react rxn1 all 1 0.0 5.0 rxn1_pre rxn1_post rxn1.map prob 0.01 1234

fix             npt_grp_react npt_grp_REACT npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0

# -- 原子轨迹转储（OVITO 兼容） --
dump            traj all custom 200 traj.lammpstrj &
                id mol type x y z vx vy vz
dump_modify     traj sort id

# -- 键拓扑转储（OVITO：作为拓扑加载） --
compute         bnd all property/local batom1 batom2 btype
dump            bonds all local 200 bonds.dump c_bnd[1] c_bnd[2] c_bnd[3]
dump_modify     bonds colname 1 batom1 colname 2 batom2 colname 3 btype

# -- 热力学输出包含来自 fix bond/react 的反应计数 --
thermo          200
thermo_style    custom step temp pe ke etotal press density f_rxns[1]
run             1000

# ====================================================================
# 阶段 5：写入最终状态
# ====================================================================
write_data      final.data
"""

from molpy.engine import LAMMPSEngine

script = mp.Script.from_text("run", lammps_script, language="other")
script.tags.add("input")

engine = LAMMPSEngine(executable="lmp_serial", check_executable=False)
proc = engine.run(
    script,
    workdir=output_dir,
    capture_output=True,
    check=False,
    timeout=600,
)

for line in proc.stdout.splitlines()[-20:]:
    print(line)
assert proc.returncode == 0, f"LAMMPS failed:\n{proc.stderr or proc.stdout[-1000:]}"

final = mp.io.read_lammps_data(output_dir / "final.data", atom_style="full")
n_atoms_final = final["atoms"].nrows
n_mols_final = len(set(final["atoms"]["mol_id"]))
print(f"final: {n_atoms_final} atoms, {n_mols_final} molecules (started with 36)")

site = find_port(left, "$")
carbon = [a for a in find_neighbors(left, site) if a.get("element") == "C"][0]
print(f"site: {carbon.get('element')} name={carbon.get('name')}")
for nb in find_neighbors(left, carbon):
    print(f"  neighbor: {nb.get('element')} name={nb.get('name')}")

print((output_dir / "rxn1.map").read_text()[:500])