REACTER Protocol¶

"How do I build a cross-linked polymer network that LAMMPS can react during MD?"

This guide builds a cross-linked polymer network by combining:

• MolPy (build monomers, assign types, generate reaction templates)
• LAMMPS fix bond/react (execute reactions during MD)

Reaction: Dehydration (ether formation)

• $\mathrm{R{-}OH + HO{-}R' \rightarrow R{-}O{-}R' + H_2O}$

What You Will Do¶

1. Build EO2 and EO3 monomers (with explicit reactive ports)
2. Generate fix bond/react templates for all pair combinations
3. Pack a mixed EO2/EO3 box with packmol as the starting configuration
4. Export a runnable LAMMPS input that loads templates and runs reactive MD

Overview: TemplateReacter Workflow¶

We use TemplateReacter to generate explicit before/after reaction templates and then feed those templates into LAMMPS fix bond/react.

The key difference from a linear polymer workflow is that EO3 creates branch points, so reactions can form a percolating network.

Setup: Imports¶

In [1]:

import numpy as np
import molpy as mp
from molpy.core.atomistic import Atomistic
from molpy.core.frame import Frame
from molpy.core.box import Box
from molpy.compute.rdkit import Generate3D
from molpy.adapter.rdkit import RDKitAdapter
from molpy.parser.smiles import parse_bigsmiles, bigsmilesir_to_monomer
from molpy.reacter import (
    select_port,
    select_c_neighbor,
    select_hydroxyl_group,
    select_hydroxyl_h_only,
    form_single_bond,
    find_port_atom,
)
from molpy.reacter.template import TemplateReacter, TemplateResult, write_template_files
from molpy.typifier.atomistic import OplsAtomisticTypifier
from molpy.io.data.lammps import LammpsDataWriter
from molpy.io.forcefield.lammps import LAMMPSForceFieldWriter
from molpy.pack import InsideBoxConstraint, Molpack
from pathlib import Path

import numpy as np import molpy as mp from molpy.core.atomistic import Atomistic from molpy.core.frame import Frame from molpy.core.box import Box from molpy.compute.rdkit import Generate3D from molpy.adapter.rdkit import RDKitAdapter from molpy.parser.smiles import parse_bigsmiles, bigsmilesir_to_monomer from molpy.reacter import ( select_port, select_c_neighbor, select_hydroxyl_group, select_hydroxyl_h_only, form_single_bond, find_port_atom, ) from molpy.reacter.template import TemplateReacter, TemplateResult, write_template_files from molpy.typifier.atomistic import OplsAtomisticTypifier from molpy.io.data.lammps import LammpsDataWriter from molpy.io.forcefield.lammps import LAMMPSForceFieldWriter from molpy.pack import InsideBoxConstraint, Molpack from pathlib import Path

Background: Reaction Templates for LAMMPS fix bond/react¶

LAMMPS fix bond/react needs explicit before/after templates.

In this notebook we generate templates for all monomer pairs (order can matter):

• EO2–EO2
• EO2–EO3
• EO3–EO3

Each template encodes reactants, products, atom mapping, and the leaving group ($H_2O$).

Checkpoint: you should see rxn*_{pre,post}.mol and rxn*.map files under the output directory.

Background: Packing an Initial Configuration¶

Before running reactive MD, we need a reasonable starting configuration: a periodic box filled with a mixed EO2/EO3 melt.

We will control:

• the EO2:EO3 composition
• the target density (sets the box size)
• Packmol overlap constraints

Checkpoint: Packmol finishes and the packed box contains the requested numbers of each monomer.

Step 1: Load Force Field and Create Typifier¶

Load a force field (OPLS-AA here) and create a typifier.

Checkpoint: after typing, atoms/bonds have assigned types suitable for LAMMPS export.

In [2]:

# Load force field
forcefield_path = "oplsaa.xml"
ff = mp.io.read_xml_forcefield(forcefield_path)
typifier = OplsAtomisticTypifier(ff, strict_typing=False)

print("✅ Force field loaded successfully")

# Load force field forcefield_path = "oplsaa.xml" ff = mp.io.read_xml_forcefield(forcefield_path) typifier = OplsAtomisticTypifier(ff, strict_typing=False) print("✅ Force field loaded successfully")

2025-12-30 16:08:58,241 - molpy.potential.dihedral.opls - WARNING - RB coefficients do not lie on the ideal 4-term OPLS manifold (C0+C1+C2+C3+C4 = 10.041600, expected ≈ 0). Conversion will preserve forces and relative energies exactly, but will introduce a constant energy offset of ΔE = 10.041600 kJ/mol. This does not affect MD simulations.

✅ Force field loaded successfully

Step 2: Define Helper Functions¶

We define small helpers to keep the notebook readable:

• build_monomer(...): build a monomer from BigSMILES and assign 3D coordinates + types
• run_reaction_with_template(...): run the dehydration reaction and capture a fix bond/react template

In [3]:

def build_monomer(bigsmiles: str, typifier: OplsAtomisticTypifier) -> Atomistic:
    """Build a monomer from BigSMILES string with 3D coordinates and types."""
    ir = parse_bigsmiles(bigsmiles)
    monomer = bigsmilesir_to_monomer(ir)

    adapter = RDKitAdapter(internal=monomer)
    generate_3d = Generate3D(add_hydrogens=True, embed=True, optimize=True, update_internal=True)
    adapter = generate_3d(adapter)
    monomer = adapter.get_internal()

    monomer.get_topo(gen_angle=True, gen_dihe=True)

    for idx, atom in enumerate(monomer.atoms):
        atom["id"] = idx + 1

    typifier.typify(monomer)

    return monomer


def run_reaction_with_template(
    left: Atomistic,
    right: Atomistic,
    port_L: str,
    port_R: str,
    reaction_name: str = "crosslinking",
    radius: int = 6,
 ) -> TemplateResult:
    """Run dehydration reaction and generate template using TemplateReacter.

    Note: If both left and right have the same port name (e.g., both have '$'),
    this function will automatically select different ports:
    - port_L will use the rightmost port (higher atom id)
    - port_R will use the leftmost port (lower atom id)
    """
    template_reacter = TemplateReacter(
        name=reaction_name,
        anchor_selector_left=select_c_neighbor,
        anchor_selector_right=select_port,
        leaving_selector_left=select_hydroxyl_group,
        leaving_selector_right=select_hydroxyl_h_only,
        bond_former=form_single_bond,
        radius=radius,
    )

    if port_L == port_R and port_L in [a.get("port") for a in left.atoms] and port_L in [a.get("port") for a in right.atoms]:
        left_ports = [a for a in left.atoms if a.get("port") == port_L]
        right_ports = [a for a in right.atoms if a.get("port") == port_R]

        if len(left_ports) >= 2 and len(right_ports) >= 1:
            port_atom_L = max(left_ports, key=lambda a: a.get("id", 0))
            port_atom_R = min(right_ports, key=lambda a: a.get("id", 0))
        elif len(left_ports) >= 1 and len(right_ports) >= 2:
            port_atom_L = min(left_ports, key=lambda a: a.get("id", 0))
            port_atom_R = max(right_ports, key=lambda a: a.get("id", 0))
        else:
            port_atom_L = find_port_atom(left, port_L)
            port_atom_R = find_port_atom(right, port_R)
    else:
        port_atom_L = find_port_atom(left, port_L)
        port_atom_R = find_port_atom(right, port_R)

    _result, template = template_reacter.run_with_template(
        left=left,
        right=right,
        port_atom_L=port_atom_L,
        port_atom_R=port_atom_R,
        compute_topology=True,
        record_intermediates=False,
    )

    return template


print("✅ Helper functions defined")

def build_monomer(bigsmiles: str, typifier: OplsAtomisticTypifier) -> Atomistic: """Build a monomer from BigSMILES string with 3D coordinates and types.""" ir = parse_bigsmiles(bigsmiles) monomer = bigsmilesir_to_monomer(ir) adapter = RDKitAdapter(internal=monomer) generate_3d = Generate3D(add_hydrogens=True, embed=True, optimize=True, update_internal=True) adapter = generate_3d(adapter) monomer = adapter.get_internal() monomer.get_topo(gen_angle=True, gen_dihe=True) for idx, atom in enumerate(monomer.atoms): atom["id"] = idx + 1 typifier.typify(monomer) return monomer def run_reaction_with_template( left: Atomistic, right: Atomistic, port_L: str, port_R: str, reaction_name: str = "crosslinking", radius: int = 6, ) -> TemplateResult: """Run dehydration reaction and generate template using TemplateReacter. Note: If both left and right have the same port name (e.g., both have '$'), this function will automatically select different ports: - port_L will use the rightmost port (higher atom id) - port_R will use the leftmost port (lower atom id) """ template_reacter = TemplateReacter( name=reaction_name, anchor_selector_left=select_c_neighbor, anchor_selector_right=select_port, leaving_selector_left=select_hydroxyl_group, leaving_selector_right=select_hydroxyl_h_only, bond_former=form_single_bond, radius=radius, ) if port_L == port_R and port_L in [a.get("port") for a in left.atoms] and port_L in [a.get("port") for a in right.atoms]: left_ports = [a for a in left.atoms if a.get("port") == port_L] right_ports = [a for a in right.atoms if a.get("port") == port_R] if len(left_ports) >= 2 and len(right_ports) >= 1: port_atom_L = max(left_ports, key=lambda a: a.get("id", 0)) port_atom_R = min(right_ports, key=lambda a: a.get("id", 0)) elif len(left_ports) >= 1 and len(right_ports) >= 2: port_atom_L = min(left_ports, key=lambda a: a.get("id", 0)) port_atom_R = max(right_ports, key=lambda a: a.get("id", 0)) else: port_atom_L = find_port_atom(left, port_L) port_atom_R = find_port_atom(right, port_R) else: port_atom_L = find_port_atom(left, port_L) port_atom_R = find_port_atom(right, port_R) _result, template = template_reacter.run_with_template( left=left, right=right, port_atom_L=port_atom_L, port_atom_R=port_atom_R, compute_topology=True, record_intermediates=False, ) return template print("✅ Helper functions defined")

✅ Helper functions defined

Template generation helper¶

run_reaction_with_template(...) wraps TemplateReacter.run_with_template(...).

Checkpoint: running it on a monomer pair yields a template that can be written to *.mol and *.map files.

Step 3: Build EO2 and EO3 Monomers¶

We build two monomers from BigSMILES and verify their reactive ports:

• EO2: linear (2 ports)
• EO3: branched (3 ports)

Checkpoint: EO2 and EO3 build successfully and their port labels are detected.

In [4]:

def build_eo2_monomer(typifier: OplsAtomisticTypifier) -> Atomistic:
    """Build EO2 monomer: HO-CH2CH2-O-CH2CH2-OH"""
    bigsmiles = "{[][$]OCCOCCOCCO[$][]}"
    return build_monomer(bigsmiles, typifier)


def build_eo3_monomer(typifier: OplsAtomisticTypifier) -> Atomistic:
    """Build EO3 monomer: 3-arm branched structure C(COCCO[$])(COCCO[$])COCCO[$]"""
    bigsmiles = "{[]C(COCCO[$])(COCCO[$])COCCO[$][]}"
    return build_monomer(bigsmiles, typifier)


# Build EO2 and EO3 monomers
eo2_monomer = build_eo2_monomer(typifier)
eo3_monomer = build_eo3_monomer(typifier)

eo2_ports = [atom.get("port") for atom in eo2_monomer.atoms if atom.get("port") is not None]
eo3_ports = [atom.get("port") for atom in eo3_monomer.atoms if atom.get("port") is not None]

print("✅ EO2 monomer built:")
print(f"   Atoms: {len(eo2_monomer.atoms)}")
print(f"   Ports: {eo2_ports}")

print("✅ EO3 monomer built:")
print(f"   Atoms: {len(eo3_monomer.atoms)}")
print(f"   Ports: {eo3_ports}")

def build_eo2_monomer(typifier: OplsAtomisticTypifier) -> Atomistic: """Build EO2 monomer: HO-CH2CH2-O-CH2CH2-OH""" bigsmiles = "{[][$]OCCOCCOCCO[$][]}" return build_monomer(bigsmiles, typifier) def build_eo3_monomer(typifier: OplsAtomisticTypifier) -> Atomistic: """Build EO3 monomer: 3-arm branched structure C(COCCO[$])(COCCO[$])COCCO[$]""" bigsmiles = "{[]C(COCCO[$])(COCCO[$])COCCO[$][]}" return build_monomer(bigsmiles, typifier) # Build EO2 and EO3 monomers eo2_monomer = build_eo2_monomer(typifier) eo3_monomer = build_eo3_monomer(typifier) eo2_ports = [atom.get("port") for atom in eo2_monomer.atoms if atom.get("port") is not None] eo3_ports = [atom.get("port") for atom in eo3_monomer.atoms if atom.get("port") is not None] print("✅ EO2 monomer built:") print(f" Atoms: {len(eo2_monomer.atoms)}") print(f" Ports: {eo2_ports}") print("✅ EO3 monomer built:") print(f" Atoms: {len(eo3_monomer.atoms)}") print(f" Ports: {eo3_ports}")

/opt/buildhome/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/compute/rdkit.py:158: UserWarning: UFF optimization returned code 1. Code 1 typically means convergence not reached within 200 iterations. The structure may still be improved.
  warnings.warn(msg, UserWarning)

✅ EO2 monomer built:
   Atoms: 24
   Ports: ['$', '$']
✅ EO3 monomer built:
   Atoms: 38
   Ports: ['$', '$', '$']

Step 4: Generate fix bond/react Templates¶

Generate templates for EO2/EO3 reactions and write them to the output directory.

Since both EO2 and EO3 have identical reactive ports ($) and the dehydration reaction is symmetric, we only need three templates:

• EO2 + EO2 (rxn1) - symmetric reaction
• EO2 + EO3 (rxn2) - covers both EO2+EO3 and EO3+EO2 directions
• EO3 + EO3 (rxn3) - symmetric reaction

Note: LAMMPS fix bond/react can match templates bidirectionally for symmetric reactions, so we don't need separate templates for EO2+EO3 and EO3+EO2.

Checkpoint: each rxn* writes a *_pre.mol, *_post.mol, and *.map file.

In [5]:

# Create output directory
output_dir = Path("user-guide-output") / "04_polymer_crosslinking"
output_dir.mkdir(parents=True, exist_ok=True)
print(f"Outputs: {output_dir.resolve()}")

# Generate reaction templates
templates: list[tuple[str, TemplateResult]] = []

# Template 1: EO2 + EO2
print("\nGenerating reaction template: EO2 + EO2...")
template1 = run_reaction_with_template(
    left=eo2_monomer.copy(),
    right=eo2_monomer.copy(),
    port_L="$",
    port_R="$",
    reaction_name="eo2_eo2",
    radius=4,
)
write_template_files(
    base_path=output_dir / "rxn1",
    template=template1,
    typifier=typifier,
)
templates.append(("rxn1", template1))
print("✅ Reaction template rxn1 generated")

# Template 2: EO2 + EO3
print("\nGenerating reaction template: EO2 + EO3...")
template2 = run_reaction_with_template(
    left=eo2_monomer.copy(),
    right=eo3_monomer.copy(),
    port_L="$",
    port_R="$",
    reaction_name="eo2_eo3",
    radius=4,
)
write_template_files(
    base_path=output_dir / "rxn2",
    template=template2,
    typifier=typifier,
)
templates.append(("rxn2", template2))
print("✅ Reaction template rxn2 generated")

# Template 3: EO3 + EO3 (symmetric reaction)
print("\nGenerating reaction template: EO3 + EO2...")
template3 = run_reaction_with_template(
    left=eo3_monomer.copy(),
    right=eo3_monomer.copy(),
    port_L="$",
    port_R="$",
    reaction_name="eo3_eo3",
    radius=4,
)
base_path = (output_dir / "rxn3",)
templates.append(("rxn3", template3))
print("✅ Reaction template rxn3 generated")

# Print summary for all templates
print("\n✅ All reaction templates generated:")
for rxn_name, template in templates:
    print(f"   - {rxn_name}_pre.mol, {rxn_name}_post.mol, {rxn_name}.map")
    print(
        f"     Pre: {len(list(template.pre.atoms))} atoms, Post: {len(list(template.post.atoms))} atoms"
    )

# Load template frames for type collection
template_frames: list[Frame] = []
for _rxn_name, template in templates:
    pre_frame = template.pre.to_frame()
    post_frame = template.post.to_frame()
    template_frames.extend([pre_frame, post_frame])

# Create output directory output_dir = Path("user-guide-output") / "04_polymer_crosslinking" output_dir.mkdir(parents=True, exist_ok=True) print(f"Outputs: {output_dir.resolve()}") # Generate reaction templates templates: list[tuple[str, TemplateResult]] = [] # Template 1: EO2 + EO2 print("\nGenerating reaction template: EO2 + EO2...") template1 = run_reaction_with_template( left=eo2_monomer.copy(), right=eo2_monomer.copy(), port_L="$", port_R="$", reaction_name="eo2_eo2", radius=4, ) write_template_files( base_path=output_dir / "rxn1", template=template1, typifier=typifier, ) templates.append(("rxn1", template1)) print("✅ Reaction template rxn1 generated") # Template 2: EO2 + EO3 print("\nGenerating reaction template: EO2 + EO3...") template2 = run_reaction_with_template( left=eo2_monomer.copy(), right=eo3_monomer.copy(), port_L="$", port_R="$", reaction_name="eo2_eo3", radius=4, ) write_template_files( base_path=output_dir / "rxn2", template=template2, typifier=typifier, ) templates.append(("rxn2", template2)) print("✅ Reaction template rxn2 generated") # Template 3: EO3 + EO3 (symmetric reaction) print("\nGenerating reaction template: EO3 + EO2...") template3 = run_reaction_with_template( left=eo3_monomer.copy(), right=eo3_monomer.copy(), port_L="$", port_R="$", reaction_name="eo3_eo3", radius=4, ) base_path = (output_dir / "rxn3",) templates.append(("rxn3", template3)) print("✅ Reaction template rxn3 generated") # Print summary for all templates print("\n✅ All reaction templates generated:") for rxn_name, template in templates: print(f" - {rxn_name}_pre.mol, {rxn_name}_post.mol, {rxn_name}.map") print( f" Pre: {len(list(template.pre.atoms))} atoms, Post: {len(list(template.post.atoms))} atoms" ) # Load template frames for type collection template_frames: list[Frame] = [] for _rxn_name, template in templates: pre_frame = template.pre.to_frame() post_frame = template.post.to_frame() template_frames.extend([pre_frame, post_frame])

Outputs: /opt/buildhome/repo/docs/user-guide/user-guide-output/04_polymer_crosslinking

Generating reaction template: EO2 + EO2...

✅ Written: rxn1_pre.mol, rxn1_post.mol, rxn1.map
✅ Reaction template rxn1 generated

Generating reaction template: EO2 + EO3...

✅ Written: rxn2_pre.mol, rxn2_post.mol, rxn2.map
✅ Reaction template rxn2 generated

Generating reaction template: EO3 + EO2...
✅ Reaction template rxn3 generated

✅ All reaction templates generated:
   - rxn1_pre.mol, rxn1_post.mol, rxn1.map
     Pre: 23 atoms, Post: 23 atoms
   - rxn2_pre.mol, rxn2_post.mol, rxn2.map
     Pre: 23 atoms, Post: 23 atoms
   - rxn3_pre.mol, rxn3_post.mol, rxn3.map
     Pre: 23 atoms, Post: 23 atoms

Step 5: Pack the Initial EO2/EO3 Mixture¶

Pack a mixed EO2/EO3 system into a periodic box using Packmol (via Molpack).

Checkpoint: packing completes and you get a LAMMPS *.data / *.ff pair under the output directory.

In [6]:

def calculate_molecular_weight(monomer: Atomistic) -> float:
    """Calculate molecular weight of a monomer."""
    from molpy.core.element import Element

    total_mw = 0.0
    for atom in monomer.atoms:
        symbol = atom.get("symbol", "C")
        symbol_upper = symbol.upper() if symbol else "C"
        element = Element(symbol_upper)
        total_mw += element.mass
    return total_mw


def calculate_box_size_from_density(
    n_monomers: int,
    monomer: Atomistic,
    target_density: float = 1.0,
 ) -> float:
    """Calculate box length (Å) from target density (g/cm^3)."""
    mw = calculate_molecular_weight(monomer)
    total_mw = n_monomers * mw
    NA = 6.022e23
    total_mass_g = total_mw / NA
    volume_cm3 = total_mass_g / target_density
    volume_angstrom3 = volume_cm3 * 1e24
    box_length_angstrom = volume_angstrom3 ** (1.0 / 3.0)
    return box_length_angstrom


def collect_types_from_frames(*frames: Frame) -> dict[str, set[str]]:
    """Collect all type names from multiple frames."""
    atom_types: set[str] = set()
    bond_types: set[str] = set()
    angle_types: set[str] = set()
    dihedral_types: set[str] = set()

    for frame in frames:
        if "atoms" in frame and "type" in frame["atoms"]:
            for atom_type in frame["atoms"]["type"]:
                if atom_type:
                    type_str = str(atom_type)
                    if not type_str.isdigit():
                        atom_types.add(type_str)

        if "bonds" in frame and "type" in frame["bonds"]:
            for bond_type in frame["bonds"]["type"]:
                if bond_type:
                    type_str = str(bond_type)
                    if not type_str.isdigit():
                        bond_types.add(type_str)

        if "angles" in frame and "type" in frame["angles"]:
            for angle_type in frame["angles"]["type"]:
                if angle_type:
                    type_str = str(angle_type)
                    if not type_str.isdigit():
                        angle_types.add(type_str)

        if "dihedrals" in frame and "type" in frame["dihedrals"]:
            for dihedral_type in frame["dihedrals"]["type"]:
                if dihedral_type:
                    type_str = str(dihedral_type)
                    if not type_str.isdigit():
                        dihedral_types.add(type_str)

    return {
        "atom_types": atom_types,
        "bond_types": bond_types,
        "angle_types": angle_types,
        "dihedral_types": dihedral_types,
    }


def clean_port_metadata(monomer: Atomistic) -> Atomistic:
    """Remove port metadata fields to avoid packing inconsistencies."""
    cleaned = monomer.copy()
    port_fields = [
        "port",
        "port_role",
        "port_bond_kind",
        "port_compat",
        "port_priority",
        "ports",
    ]
    for atom in cleaned.atoms:
        for field in port_fields:
            if field in atom.data:
                del atom.data[field]
    return cleaned

def calculate_molecular_weight(monomer: Atomistic) -> float: """Calculate molecular weight of a monomer.""" from molpy.core.element import Element total_mw = 0.0 for atom in monomer.atoms: symbol = atom.get("symbol", "C") symbol_upper = symbol.upper() if symbol else "C" element = Element(symbol_upper) total_mw += element.mass return total_mw def calculate_box_size_from_density( n_monomers: int, monomer: Atomistic, target_density: float = 1.0, ) -> float: """Calculate box length (Å) from target density (g/cm^3).""" mw = calculate_molecular_weight(monomer) total_mw = n_monomers * mw NA = 6.022e23 total_mass_g = total_mw / NA volume_cm3 = total_mass_g / target_density volume_angstrom3 = volume_cm3 * 1e24 box_length_angstrom = volume_angstrom3 ** (1.0 / 3.0) return box_length_angstrom def collect_types_from_frames(*frames: Frame) -> dict[str, set[str]]: """Collect all type names from multiple frames.""" atom_types: set[str] = set() bond_types: set[str] = set() angle_types: set[str] = set() dihedral_types: set[str] = set() for frame in frames: if "atoms" in frame and "type" in frame["atoms"]: for atom_type in frame["atoms"]["type"]: if atom_type: type_str = str(atom_type) if not type_str.isdigit(): atom_types.add(type_str) if "bonds" in frame and "type" in frame["bonds"]: for bond_type in frame["bonds"]["type"]: if bond_type: type_str = str(bond_type) if not type_str.isdigit(): bond_types.add(type_str) if "angles" in frame and "type" in frame["angles"]: for angle_type in frame["angles"]["type"]: if angle_type: type_str = str(angle_type) if not type_str.isdigit(): angle_types.add(type_str) if "dihedrals" in frame and "type" in frame["dihedrals"]: for dihedral_type in frame["dihedrals"]["type"]: if dihedral_type: type_str = str(dihedral_type) if not type_str.isdigit(): dihedral_types.add(type_str) return { "atom_types": atom_types, "bond_types": bond_types, "angle_types": angle_types, "dihedral_types": dihedral_types, } def clean_port_metadata(monomer: Atomistic) -> Atomistic: """Remove port metadata fields to avoid packing inconsistencies.""" cleaned = monomer.copy() port_fields = [ "port", "port_role", "port_bond_kind", "port_compat", "port_priority", "ports", ] for atom in cleaned.atoms: for field in port_fields: if field in atom.data: del atom.data[field] return cleaned

Choosing a target density¶

The target density determines the box size. For polymer melts, values around 0.8–1.2 g/cm³ are common.

If you change density (or composition), re-check the resulting box length and whether packing remains stable.

In [7]:

# Pack monomers (mix of EO2 and EO3)
n_eo2 = 27
n_eo3 = 9
n_monomers = n_eo2 + n_eo3
target_density = 1.2

# Compute box length from target density
mw_eo2 = calculate_molecular_weight(eo2_monomer)
mw_eo3 = calculate_molecular_weight(eo3_monomer)
total_mw = n_eo2 * mw_eo2 + n_eo3 * mw_eo3

NA = 6.022e23
total_mass_g = total_mw / NA
volume_cm3 = total_mass_g / target_density
volume_angstrom3 = volume_cm3 * 1e24
box_length = volume_angstrom3 ** (1.0 / 3.0)

print(f"\n✅ Calculated box size: {box_length:.2f} Å")
print(f"   Mix: {n_eo2} EO2 + {n_eo3} EO3 = {n_monomers} total monomers")

# Prepare frames for packing (clean port metadata first)
eo2_clean = clean_port_metadata(eo2_monomer)
eo3_clean = clean_port_metadata(eo3_monomer)
eo2_frame = eo2_clean.to_frame()
eo3_frame = eo3_clean.to_frame()

# Normalize charge field name for pack/export
if "charge" in eo2_frame["atoms"]:
    eo2_frame["atoms"]["q"] = eo2_frame["atoms"]["charge"]
if "charge" in eo3_frame["atoms"]:
    eo3_frame["atoms"]["q"] = eo3_frame["atoms"]["charge"]

# Use packmol for packing (Packmol-only)
pack_workdir = output_dir / "packmol_work"
packer = Molpack(workdir=pack_workdir)

origin = np.array([0.0, 0.0, 0.0])
length = np.array([box_length, box_length, box_length])
box_constraint = InsideBoxConstraint(length=length, origin=origin)

_ = packer.add_target(eo2_frame, number=n_eo2, constraint=box_constraint)
_ = packer.add_target(eo3_frame, number=n_eo3, constraint=box_constraint)

print(f"   Packing {n_eo2} EO2 + {n_eo3} EO3 monomers...")
packed_frame = packer.optimize(max_steps=10000, seed=42)
print(f"✅ Packed system: {packed_frame['atoms'].nrows} atoms")

# Ensure required charge fields exist
if "charge" in packed_frame["atoms"] and "q" not in packed_frame["atoms"]:
    packed_frame["atoms"]["q"] = packed_frame["atoms"]["charge"]
elif "q" in packed_frame["atoms"] and "charge" not in packed_frame["atoms"]:
    packed_frame["atoms"]["charge"] = packed_frame["atoms"]["q"]

# Add box metadata for downstream writers
box = Box.cubic(length=box_length)
packed_frame.metadata["box"] = box

# Collect types from config and templates
all_frames = [packed_frame] + template_frames
types_from_config = collect_types_from_frames(*all_frames)

packed_frame.metadata["type_labels"] = {
    "atom_types": sorted(list(types_from_config["atom_types"])),
    "bond_types": sorted(list(types_from_config["bond_types"])),
    "angle_types": sorted(list(types_from_config["angle_types"])),
    "dihedral_types": sorted(list(types_from_config["dihedral_types"])),
}

# Write force field and data files
ff_path = output_dir / "case4_crosslinking.ff"
ff_writer = LAMMPSForceFieldWriter(ff_path)
ff_writer.write(
    ff,
    atom_types=types_from_config["atom_types"],
    bond_types=types_from_config["bond_types"],
    angle_types=types_from_config["angle_types"],
    dihedral_types=types_from_config["dihedral_types"],
 )

data_path = output_dir / "case4_crosslinking.data"
data_writer = LammpsDataWriter(data_path, atom_style="full")
data_writer.write(packed_frame)

print(f"✅ Written: {ff_path.name}, {data_path.name}")

# Pack monomers (mix of EO2 and EO3) n_eo2 = 27 n_eo3 = 9 n_monomers = n_eo2 + n_eo3 target_density = 1.2 # Compute box length from target density mw_eo2 = calculate_molecular_weight(eo2_monomer) mw_eo3 = calculate_molecular_weight(eo3_monomer) total_mw = n_eo2 * mw_eo2 + n_eo3 * mw_eo3 NA = 6.022e23 total_mass_g = total_mw / NA volume_cm3 = total_mass_g / target_density volume_angstrom3 = volume_cm3 * 1e24 box_length = volume_angstrom3 ** (1.0 / 3.0) print(f"\n✅ Calculated box size: {box_length:.2f} Å") print(f" Mix: {n_eo2} EO2 + {n_eo3} EO3 = {n_monomers} total monomers") # Prepare frames for packing (clean port metadata first) eo2_clean = clean_port_metadata(eo2_monomer) eo3_clean = clean_port_metadata(eo3_monomer) eo2_frame = eo2_clean.to_frame() eo3_frame = eo3_clean.to_frame() # Normalize charge field name for pack/export if "charge" in eo2_frame["atoms"]: eo2_frame["atoms"]["q"] = eo2_frame["atoms"]["charge"] if "charge" in eo3_frame["atoms"]: eo3_frame["atoms"]["q"] = eo3_frame["atoms"]["charge"] # Use packmol for packing (Packmol-only) pack_workdir = output_dir / "packmol_work" packer = Molpack(workdir=pack_workdir) origin = np.array([0.0, 0.0, 0.0]) length = np.array([box_length, box_length, box_length]) box_constraint = InsideBoxConstraint(length=length, origin=origin) _ = packer.add_target(eo2_frame, number=n_eo2, constraint=box_constraint) _ = packer.add_target(eo3_frame, number=n_eo3, constraint=box_constraint) print(f" Packing {n_eo2} EO2 + {n_eo3} EO3 monomers...") packed_frame = packer.optimize(max_steps=10000, seed=42) print(f"✅ Packed system: {packed_frame['atoms'].nrows} atoms") # Ensure required charge fields exist if "charge" in packed_frame["atoms"] and "q" not in packed_frame["atoms"]: packed_frame["atoms"]["q"] = packed_frame["atoms"]["charge"] elif "q" in packed_frame["atoms"] and "charge" not in packed_frame["atoms"]: packed_frame["atoms"]["charge"] = packed_frame["atoms"]["q"] # Add box metadata for downstream writers box = Box.cubic(length=box_length) packed_frame.metadata["box"] = box # Collect types from config and templates all_frames = [packed_frame] + template_frames types_from_config = collect_types_from_frames(*all_frames) packed_frame.metadata["type_labels"] = { "atom_types": sorted(list(types_from_config["atom_types"])), "bond_types": sorted(list(types_from_config["bond_types"])), "angle_types": sorted(list(types_from_config["angle_types"])), "dihedral_types": sorted(list(types_from_config["dihedral_types"])), } # Write force field and data files ff_path = output_dir / "case4_crosslinking.ff" ff_writer = LAMMPSForceFieldWriter(ff_path) ff_writer.write( ff, atom_types=types_from_config["atom_types"], bond_types=types_from_config["bond_types"], angle_types=types_from_config["angle_types"], dihedral_types=types_from_config["dihedral_types"], ) data_path = output_dir / "case4_crosslinking.data" data_writer = LammpsDataWriter(data_path, atom_style="full") data_writer.write(packed_frame) print(f"✅ Written: {ff_path.name}, {data_path.name}")

✅ Calculated box size: 20.47 Å
   Mix: 27 EO2 + 9 EO3 = 36 total monomers
   Packing 27 EO2 + 9 EO3 monomers...

✅ Packed system: 990 atoms
✅ Written: case4_crosslinking.ff, case4_crosslinking.data

Step 6: Generate a LAMMPS Input Script¶

Generate a runnable LAMMPS input that:

• reads the packed *.data and *.ff
• loads the fix bond/react templates
• runs equilibration + reactive MD

Checkpoint: the script is written under the output directory.

In [8]:

# Build molecule and reaction commands for fix bond/react
molecule_commands = []
react_commands = []

for rxn_name, _ in templates:
    molecule_commands.append(f"molecule   {rxn_name}_pre {rxn_name}_pre.mol")
    molecule_commands.append(f"molecule   {rxn_name}_post {rxn_name}_post.mol")
    react_commands.append(
        f"  react {rxn_name} all 1 0.0 5 {rxn_name}_pre {rxn_name}_post {rxn_name}.map prob 0.2 5123 rescale_charges yes &"
    )

script_content = f"""# LAMMPS input script for Case 3: Cross-linked EO2/EO3 System
# Force field: OPLS-AA
# Reaction: Dehydration to form ether bonds
# Monomers: {n_eo2} EO2 + {n_eo3} EO3
#
units           real
atom_style      full
boundary        p p p
dimension       3

read_data       case4_crosslinking.data &
    extra/bond/per/atom 25 &
    extra/angle/per/atom 25 &
    extra/dihedral/per/atom 25 &
    extra/improper/per/atom 25 &
    extra/special/per/atom 25
include         case4_crosslinking.ff
kspace_style    pppm 1.0e-5
special_bonds   lj/coul 0.0 0.0 0.5

neighbor        2.0 bin
neigh_modify    delay 0 every 1 check yes

print "=========================================="
print "Step 1: Energy Minimization"
print "=========================================="
minimize        1.0e-4 1.0e-4 1000 10000

print "=========================================="
print "Step 2: NPT Equilibration"
print "=========================================="
velocity all create 300.0 1234
timestep        0.1
fix             shake all shake 0.0001 20 0 t opls_140 opls_155 opls_185
fix             npt all npt temp 300.0 300.0 100.0 iso 1.5 1.5 1000.0
dump            2 all custom 500 equil.dump id type x y z
run             1000
unfix           npt
undump          2

print "=========================================="
print "Step 3: Reactive MD with fix bond/react"
print "=========================================="

# Read reaction templates
{chr(10).join(molecule_commands)}

fix rxns all bond/react stabilization yes npt_grp .03 &
{chr(10).join(react_commands)}
    # End of reactions
unfix shake

fix             npt_grp_react all npt temp 300.0 300.0 100.0 iso 1.5 1.5 1000.0
thermo          100
thermo_style    custom elapsed temp pe ke etotal press vol density f_rxns[*]
thermo_modify   flush yes
compute 1 all property/local batom1 batom2 btype
dump            3 all custom 10 react.dump id mol type xu yu zu
dump            4 all local 10 topo.dump index c_1[1] c_1[2] c_1[3]
run             5000
undump         3
undump         4
unfix          rxns
unfix          npt_grp_react

print "=========================================="
print "Step 4: Production MD"
print "=========================================="
thermo_style    custom elapsed temp pe ke etotal press vol density

write_data      final_crosslinked.data
write_restart   final_crosslinked.restart

print "=========================================="
print "Simulation completed!"
print "=========================================="
"""

script_path = output_dir / "run_case4_crosslinking.in"
with script_path.open("w") as f:
    f.write(script_content)

print(f"✅ Generated LAMMPS input script: {script_path.name}")

# Build molecule and reaction commands for fix bond/react molecule_commands = [] react_commands = [] for rxn_name, _ in templates: molecule_commands.append(f"molecule {rxn_name}_pre {rxn_name}_pre.mol") molecule_commands.append(f"molecule {rxn_name}_post {rxn_name}_post.mol") react_commands.append( f" react {rxn_name} all 1 0.0 5 {rxn_name}_pre {rxn_name}_post {rxn_name}.map prob 0.2 5123 rescale_charges yes &" ) script_content = f"""# LAMMPS input script for Case 3: Cross-linked EO2/EO3 System # Force field: OPLS-AA # Reaction: Dehydration to form ether bonds # Monomers: {n_eo2} EO2 + {n_eo3} EO3 # units real atom_style full boundary p p p dimension 3 read_data case4_crosslinking.data & extra/bond/per/atom 25 & extra/angle/per/atom 25 & extra/dihedral/per/atom 25 & extra/improper/per/atom 25 & extra/special/per/atom 25 include case4_crosslinking.ff kspace_style pppm 1.0e-5 special_bonds lj/coul 0.0 0.0 0.5 neighbor 2.0 bin neigh_modify delay 0 every 1 check yes print "==========================================" print "Step 1: Energy Minimization" print "==========================================" minimize 1.0e-4 1.0e-4 1000 10000 print "==========================================" print "Step 2: NPT Equilibration" print "==========================================" velocity all create 300.0 1234 timestep 0.1 fix shake all shake 0.0001 20 0 t opls_140 opls_155 opls_185 fix npt all npt temp 300.0 300.0 100.0 iso 1.5 1.5 1000.0 dump 2 all custom 500 equil.dump id type x y z run 1000 unfix npt undump 2 print "==========================================" print "Step 3: Reactive MD with fix bond/react" print "==========================================" # Read reaction templates {chr(10).join(molecule_commands)} fix rxns all bond/react stabilization yes npt_grp .03 & {chr(10).join(react_commands)} # End of reactions unfix shake fix npt_grp_react all npt temp 300.0 300.0 100.0 iso 1.5 1.5 1000.0 thermo 100 thermo_style custom elapsed temp pe ke etotal press vol density f_rxns[*] thermo_modify flush yes compute 1 all property/local batom1 batom2 btype dump 3 all custom 10 react.dump id mol type xu yu zu dump 4 all local 10 topo.dump index c_1[1] c_1[2] c_1[3] run 5000 undump 3 undump 4 unfix rxns unfix npt_grp_react print "==========================================" print "Step 4: Production MD" print "==========================================" thermo_style custom elapsed temp pe ke etotal press vol density write_data final_crosslinked.data write_restart final_crosslinked.restart print "==========================================" print "Simulation completed!" print "==========================================" """ script_path = output_dir / "run_case4_crosslinking.in" with script_path.open("w") as f: f.write(script_content) print(f"✅ Generated LAMMPS input script: {script_path.name}")

✅ Generated LAMMPS input script: run_case4_crosslinking.in

Summary and Next Steps¶

You now have a complete workflow for building a cross-linked EO2/EO3 system for reactive MD:

1. Load OPLS-AA and set up typing
2. Build EO2 and EO3 monomers with reactive ports
3. Generate fix bond/react templates with TemplateReacter
4. Pack an initial mixed configuration with Packmol
5. Export LAMMPS-ready inputs and a runnable script

Generated files (under the output directory):

• LAMMPS input script: run_case4_crosslinking.in
• Reaction templates: rxn1–rxn3 (*_pre.mol, *_post.mol, *.map)
• Initial configuration: case4_crosslinking.data, case4_crosslinking.ff