Crosslinked Networks¶

This guide develops the inputs required for reactive polymer simulations with LAMMPS fix bond/react: local reaction templates, packed initial configurations, and exported topology and force-field files.

!!! note "Prerequisites" This guide requires RDKit, Packmol, and the oplsaa.xml force field. Familiarity with Stepwise Polymer Construction is assumed.

Reactive MD requires four outputs¶

LAMMPS fix bond/react drives bond formation during MD. It needs reaction templates that capture the local topology before and after each bond formation event, a packed configuration at realistic density, force field coefficients that cover every type in both the initial configuration and the post-reaction templates, and an input script that sequences these ingredients through minimisation, equilibration, and reactive MD. MolPy's BondReactReacter generates the template files, Molpack handles packing, and write_lammps_bond_react_system exports the data, force field, and templates with unified type numbering in one call.

Two monomers share one reaction template¶

We use two monomers: EO2 (linear, two $ ports) and EO3 (branched, three $ ports). Although the system contains two different species, it only needs one reaction template.

The reason is that fix bond/react matches templates by local topology, not by whole-molecule identity. The template only captures a few bond-hops around the reaction site (controlled by the radius parameter). Both EO2 and EO3 share the same arm chemistry — every reactive port sits at the end of a ...COCCO[$] segment. When the BFS stops at radius=4, it has not yet reached the EO3 branching carbon, so the extracted subgraph is structurally identical regardless of whether the port belongs to an EO2 or an EO3 molecule.

EO2:   HO–OCCOCCOCCO–OH           ← 2 identical arms
            [$]          [$]

EO3:        COCCO–OH               ← 3 identical arms
           /     [$]
     C–COCCO–OH
           \     [$]
            COCCO–OH
                  [$]

Template radius captures only the port neighbourhood:
     ...COCCO–OH   (same for every arm)
         ^^^^
         radius=4 stops here

If you had monomers with chemically different port environments (e.g. an amine reacting with an epoxide), you would need separate templates — one per distinct reaction type.

In [1]:

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 OplsAtomisticTypifier

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

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 OplsAtomisticTypifier ff = mp.io.read_xml_forcefield("oplsaa.xml") typifier = OplsAtomisticTypifier(ff, strict_typing=False)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[1], line 3
      1 from pathlib import Path
      2 import numpy as np
----> 3 import molpy as mp
      4 from molpy.core.atomistic import Atom, Atomistic
      5 from molpy.core.element import Element
      6 from molpy.typifier import OplsAtomisticTypifier

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/__init__.py:4
      1 """MolPy — Composable molecular modeling in Python."""
      3 # Submodules
----> 4 from . import adapter, data, engine, io, legacy, parser, potential, typifier
      6 # Core
      7 from .core.atomistic import Angle, Atom, Atomistic, Bond, Dihedral, Improper

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/adapter/__init__.py:12
     10 # Optional RDKit adapter
     11 try:  # pragma: no cover
---> 12     from .rdkit import MP_ID, Generate3D, OptimizeGeometry, RDKitAdapter, generate_3d
     14     _HAS_RDKIT = True
     15 except ModuleNotFoundError:  # rdkit missing

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/adapter/rdkit.py:18
     15 from rdkit import Chem
     16 from rdkit.Chem import AllChem
---> 18 from molpy.core.atomistic import Atomistic
     20 from .base import Adapter
     22 MP_ID = "mp_id"

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/core/__init__.py:1
----> 1 from .atomistic import Angle, Atom, Atomistic, Bond, Dihedral, Improper
      2 from .box import Box
      3 from .cg import Bead, CGBond, CoarseGrain

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/core/atomistic.py:232
    228         """Third substituent atom."""
    229         return self.endpoints[3]
--> 232 class Atomistic(Struct, MembershipMixin, SpatialMixin, ConnectivityMixin, molrs.Graph):
    233     """All-atom molecular structure with full topological information.
    234 
    235     Atomistic is the primary container for molecular systems in MolPy. It
   (...)    245         Atom, Bond, Angle, Dihedral, Struct, Frame
    246     """
    248     def __init__(self, **props) -> None:

AttributeError: module 'molrs' has no attribute 'Graph'

Each monomer is parsed from BigSMILES, expanded to 3D with hydrogens, and typified with OPLS-AA. Atom IDs must be assigned before typification because the force field writer expects them.

In [2]:

eo2 = mp.parser.parse_monomer("{[][$]OCCOCCOCCO[$][]}")
eo2 = mp.adapter.generate_3d(eo2, 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)

eo2 = mp.parser.parse_monomer("{[][$]OCCOCCOCCO[$][]}") eo2 = mp.adapter.generate_3d(eo2, 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)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 1
----> 1 eo2 = mp.parser.parse_monomer("{[][$]OCCOCCOCCO[$][]}")
      2 eo2 = mp.adapter.generate_3d(eo2, add_hydrogens=True, optimize=True)
      3 eo2 = eo2.get_topo(gen_angle=True, gen_dihe=True)
      4 for idx, atom in enumerate(eo2.atoms, start=1):

NameError: name 'mp' is not defined

In [3]:

eo3 = mp.parser.parse_monomer("{[]C(COCCO[$])(COCCO[$])COCCO[$][]}")
eo3 = mp.adapter.generate_3d(eo3, 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)

eo3 = mp.parser.parse_monomer("{[]C(COCCO[$])(COCCO[$])COCCO[$][]}") eo3 = mp.adapter.generate_3d(eo3, 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)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 1
----> 1 eo3 = mp.parser.parse_monomer("{[]C(COCCO[$])(COCCO[$])COCCO[$][]}")
      2 eo3 = mp.adapter.generate_3d(eo3, add_hydrogens=True, optimize=True)
      3 eo3 = eo3.get_topo(gen_angle=True, gen_dihe=True)
      4 for idx, atom in enumerate(eo3.atoms, start=1):

NameError: name 'mp' is not defined

Verify that the port markers survived 3D generation:

In [4]:

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}")

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}")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 1
----> 1 for label, mon in [("EO2", eo2), ("EO3", eo3)]:
      2     ports = [a.get("port") for a in mon.atoms if a.get("port")]
      3     print(f"{label}: atoms={len(mon.atoms)}, ports={ports}")

NameError: name 'eo2' is not defined

Template generation captures the local reaction environment¶

BondReactReacter extends Reacter with subgraph extraction. It runs a representative reaction between two monomers, then extracts the local environment around each reaction site (controlled by radius) to produce pre-reaction and post-reaction molecule templates plus an atom equivalence map.

The radius parameter controls how many bond-hops the BFS extends from the anchor atoms. A value of 4 gives enough buffer so that all type-changed atoms are at least two bonds from the template edge — a requirement of LAMMPS fix bond/react.

The reaction follows dehydration condensation, the same chemistry as in Stepwise Polymer Construction. On the left side, the anchor is the carbon neighbor of the port atom, and the leaving group is the hydroxyl (OH). On the right side, the anchor is the port atom itself, and the leaving group is one hydrogen.

In [5]:

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


def select_hydroxyl_group(struct: Atomistic, site: Atom) -> list[Atom]:
    """Find and return [O, H] — the hydroxyl leaving group."""
    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")

from molpy.reacter import ( BondReactReacter, find_neighbors, find_port, form_single_bond, select_hydrogens, select_neighbor, select_self, ) def select_hydroxyl_group(struct: Atomistic, site: Atom) -> list[Atom]: """Find and return [O, H] — the hydroxyl leaving group.""" 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")

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[5], line 1
----> 1 from molpy.reacter import (
      2     BondReactReacter,
      3     find_neighbors,
      4     find_port,

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/__init__.py:4
      1 """MolPy — Composable molecular modeling in Python."""
      3 # Submodules
----> 4 from . import adapter, data, engine, io, legacy, parser, potential, typifier
      6 # Core
      7 from .core.atomistic import Angle, Atom, Atomistic, Bond, Dihedral, Improper

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/adapter/__init__.py:12
     10 # Optional RDKit adapter
     11 try:  # pragma: no cover
---> 12     from .rdkit import MP_ID, Generate3D, OptimizeGeometry, RDKitAdapter, generate_3d
     14     _HAS_RDKIT = True
     15 except ModuleNotFoundError:  # rdkit missing

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/adapter/rdkit.py:18
     15 from rdkit import Chem
     16 from rdkit.Chem import AllChem
---> 18 from molpy.core.atomistic import Atomistic
     20 from .base import Adapter
     22 MP_ID = "mp_id"

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/core/__init__.py:1
----> 1 from .atomistic import Angle, Atom, Atomistic, Bond, Dihedral, Improper
      2 from .box import Box
      3 from .cg import Bead, CGBond, CoarseGrain

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/core/atomistic.py:232
    228         """Third substituent atom."""
    229         return self.endpoints[3]
--> 232 class Atomistic(Struct, MembershipMixin, SpatialMixin, ConnectivityMixin, molrs.Graph):
    233     """All-atom molecular structure with full topological information.
    234 
    235     Atomistic is the primary container for molecular systems in MolPy. It
   (...)    245         Atom, Bond, Angle, Dihedral, Struct, Frame
    246     """
    248     def __init__(self, **props) -> None:

AttributeError: module 'molrs' has no attribute 'Graph'

The typifier must be passed to run(). Without it, the newly formed C–O bond has no force field type and gets silently dropped from the post template — LAMMPS would then remove leaving-group atoms but never create the crosslink.

In [6]:

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

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 reacter = BondReactReacter(
      2     name="rxn1",
      3     anchor_selector_left=select_neighbor("C"),
      4     anchor_selector_right=select_self,

NameError: name 'BondReactReacter' is not defined

The pre and post templates have the same number of atoms. Atoms in the leaving group are marked for deletion in the .map file — LAMMPS removes them after applying the post-template topology.

In [7]:

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

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 1
----> 1 print(f"pre:  {len(template.pre.atoms)} atoms")
      2 print(f"post: {len(template.post.atoms)} atoms")

NameError: name 'template' is not defined

Packing at realistic density¶

Compute the box size from total molecular weight and target density, then let Packmol find non-overlapping placements. 27 bifunctional + 9 trifunctional monomers give up to 81 reactive ports.

In [8]:

from molpy.pack import InsideBoxConstraint, Molpack

N_EO2, N_EO3 = 27, 9
TARGET_DENSITY = 1.1  # g/cm³ (amorphous 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)

from molpy.pack import InsideBoxConstraint, Molpack N_EO2, N_EO3 = 27, 9 TARGET_DENSITY = 1.1 # g/cm³ (amorphous 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)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[8], line 1
----> 1 from molpy.pack import InsideBoxConstraint, Molpack
      2 
      3 N_EO2, N_EO3 = 27, 9
      4 TARGET_DENSITY = 1.1  # g/cm³ (amorphous PEO ≈ 1.1–1.2)

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/__init__.py:4
      1 """MolPy — Composable molecular modeling in Python."""
      3 # Submodules
----> 4 from . import adapter, data, engine, io, legacy, parser, potential, typifier
      6 # Core
      7 from .core.atomistic import Angle, Atom, Atomistic, Bond, Dihedral, Improper

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/adapter/__init__.py:12
     10 # Optional RDKit adapter
     11 try:  # pragma: no cover
---> 12     from .rdkit import MP_ID, Generate3D, OptimizeGeometry, RDKitAdapter, generate_3d
     14     _HAS_RDKIT = True
     15 except ModuleNotFoundError:  # rdkit missing

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/adapter/rdkit.py:18
     15 from rdkit import Chem
     16 from rdkit.Chem import AllChem
---> 18 from molpy.core.atomistic import Atomistic
     20 from .base import Adapter
     22 MP_ID = "mp_id"

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/core/__init__.py:1
----> 1 from .atomistic import Angle, Atom, Atomistic, Bond, Dihedral, Improper
      2 from .box import Box
      3 from .cg import Bead, CGBond, CoarseGrain

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/core/atomistic.py:232
    228         """Third substituent atom."""
    229         return self.endpoints[3]
--> 232 class Atomistic(Struct, MembershipMixin, SpatialMixin, ConnectivityMixin, molrs.Graph):
    233     """All-atom molecular structure with full topological information.
    234 
    235     Atomistic is the primary container for molecular systems in MolPy. It
   (...)    245         Atom, Bond, Angle, Dihedral, Struct, Frame
    246     """
    248     def __init__(self, **props) -> None:

AttributeError: module 'molrs' has no attribute 'Graph'

In [9]:

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

packed = packer.optimize(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")

packer = Molpack(workdir=Path("04_output/packmol")) constraint = InsideBoxConstraint( length=np.array([box_length] * 3), origin=np.zeros(3), ) packer.add_target(eo2.to_frame(), number=N_EO2, constraint=constraint) packer.add_target(eo3.to_frame(), number=N_EO3, constraint=constraint) packed = packer.optimize(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")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[9], line 1
----> 1 packer = Molpack(workdir=Path("04_output/packmol"))
      2 constraint = InsideBoxConstraint(
      3     length=np.array([box_length] * 3),
      4     origin=np.zeros(3),

NameError: name 'Molpack' is not defined

One call exports data, force field, and templates together¶

write_lammps_bond_react_system collects all atom, bond, angle, and dihedral types across the packed frame and every template, builds a unified type map, and writes everything to one directory. This guarantees that the numeric type IDs in 04.data, 04.ff, rxn1_pre.mol, and rxn1_post.mol are mutually consistent.

In [10]:

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},
)

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}, )

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[10], line 6
      2 
      3 output_dir = Path("04_output")
      4 output_dir.mkdir(exist_ok=True)
      5 
----> 6 mp.io.write_lammps_bond_react_system(
      7     output_dir,
      8     packed,
      9     ff,

NameError: name 'mp' is not defined

The directory now contains the packed system configuration (04.data), force field coefficients covering every type that appears in the initial state or either template (04.ff), the pre-reaction and post-reaction molecule templates (rxn1_pre.mol and rxn1_post.mol), and the atom equivalence, edge, and delete ID map (rxn1.map).

The LAMMPS input script runs a five-stage protocol¶

The simulation progresses through five stages in sequence. Stage 1 removes steric clashes from the packed configuration with conjugate-gradient energy minimisation. Stage 2 heats the system to 300 K under NVT for 5 ps so that kinetic energy is distributed before the box is allowed to relax. Stage 3 switches to NPT at 1 atm for another 5 ps to bring the density to its equilibrium value. Stage 4 activates fix bond/react: the stabilization yes keyword places freshly reacted atoms into a separate thermostat group for a brief settling period, and molecule inter restricts reactions to atoms on different molecules, preventing intramolecular ring closure. Stage 5 cleans up, recalculates molecule IDs from the new bond topology, and writes the final snapshot.

In [11]:

lammps_script = """\
# ====================================================================
# PEO crosslinked network – OPLS-AA / fix bond/react
# ====================================================================
units           real
atom_style      full
boundary        p p p

read_data       04.data
include         04.ff

# -- Long-range electrostatics --
kspace_style    pppm 1.0e-4

# -- OPLS-AA 1-4 scaling --
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

# ====================================================================
# Stage 1: Energy minimisation (remove packing overlaps)
# ====================================================================
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

# ====================================================================
# Stage 2: NVT equilibration (300 K, 5 ps)
# ====================================================================
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             5000
unfix           nvt_eq

# ====================================================================
# Stage 3: NPT equilibration (300 K, 1 atm, 5 ps)
# ====================================================================
fix             equil all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0
run             200000
run             200000
run             100000
unfix           equil

# ====================================================================
# Stage 4: Reactive MD with fix bond/react (300 K NPT, 10 ps)
# ====================================================================
molecule        rxn1_pre rxn1_pre.mol
molecule        rxn1_post rxn1_post.mol

# bond/react: attempt every step, 5 Å cutoff, intermolecular only
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

# -- Atom trajectory dump (OVITO-compatible) --
dump            traj all custom 200 traj.lammpstrj &
                id mol type x y z vx vy vz
dump_modify     traj sort id

# -- Bond topology dump (OVITO: load as topology) --
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

# -- Thermo includes reaction count from fix bond/react --
thermo          200
thermo_style    custom step temp pe ke etotal press density f_rxns[1]
run             500000

# ====================================================================
# Stage 5: Write final state
# ====================================================================
write_data      final.data
"""

lammps_script = """\ # ==================================================================== # PEO crosslinked network – OPLS-AA / fix bond/react # ==================================================================== units real atom_style full boundary p p p read_data 04.data include 04.ff # -- Long-range electrostatics -- kspace_style pppm 1.0e-4 # -- OPLS-AA 1-4 scaling -- 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 # ==================================================================== # Stage 1: Energy minimisation (remove packing overlaps) # ==================================================================== 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 # ==================================================================== # Stage 2: NVT equilibration (300 K, 5 ps) # ==================================================================== 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 5000 unfix nvt_eq # ==================================================================== # Stage 3: NPT equilibration (300 K, 1 atm, 5 ps) # ==================================================================== fix equil all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 run 200000 run 200000 run 100000 unfix equil # ==================================================================== # Stage 4: Reactive MD with fix bond/react (300 K NPT, 10 ps) # ==================================================================== molecule rxn1_pre rxn1_pre.mol molecule rxn1_post rxn1_post.mol # bond/react: attempt every step, 5 Å cutoff, intermolecular only 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 # -- Atom trajectory dump (OVITO-compatible) -- dump traj all custom 200 traj.lammpstrj & id mol type x y z vx vy vz dump_modify traj sort id # -- Bond topology dump (OVITO: load as topology) -- 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 # -- Thermo includes reaction count from fix bond/react -- thermo 200 thermo_style custom step temp pe ke etotal press density f_rxns[1] run 500000 # ==================================================================== # Stage 5: Write final state # ==================================================================== write_data final.data """

The reset_mol_ids all command at the end recalculates molecule IDs from the current bond topology. If crosslinking succeeded, the original 36 molecules collapse into a small number of connected components.

LAMMPSEngine wraps the script into a Script object and runs lmp_serial in the output directory where all data and template files already sit.

In [12]:

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:]}"

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:]}"

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[12], line 1
----> 1 from molpy.engine import LAMMPSEngine
      2 
      3 script = mp.Script.from_text("run", lammps_script, language="other")
      4 script.tags.add("input")

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/__init__.py:4
      1 """MolPy — Composable molecular modeling in Python."""
      3 # Submodules
----> 4 from . import adapter, data, engine, io, legacy, parser, potential, typifier
      6 # Core
      7 from .core.atomistic import Angle, Atom, Atomistic, Bond, Dihedral, Improper

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/adapter/__init__.py:12
     10 # Optional RDKit adapter
     11 try:  # pragma: no cover
---> 12     from .rdkit import MP_ID, Generate3D, OptimizeGeometry, RDKitAdapter, generate_3d
     14     _HAS_RDKIT = True
     15 except ModuleNotFoundError:  # rdkit missing

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/adapter/rdkit.py:18
     15 from rdkit import Chem
     16 from rdkit.Chem import AllChem
---> 18 from molpy.core.atomistic import Atomistic
     20 from .base import Adapter
     22 MP_ID = "mp_id"

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/core/__init__.py:1
----> 1 from .atomistic import Angle, Atom, Atomistic, Bond, Dihedral, Improper
      2 from .box import Box
      3 from .cg import Bead, CGBond, CoarseGrain

File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/core/atomistic.py:232
    228         """Third substituent atom."""
    229         return self.endpoints[3]
--> 232 class Atomistic(Struct, MembershipMixin, SpatialMixin, ConnectivityMixin, molrs.Graph):
    233     """All-atom molecular structure with full topological information.
    234 
    235     Atomistic is the primary container for molecular systems in MolPy. It
   (...)    245         Atom, Bond, Angle, Dihedral, Struct, Frame
    246     """
    248     def __init__(self, **props) -> None:

AttributeError: module 'molrs' has no attribute 'Graph'

!!! tip "Visualising in OVITO" Load traj.lammpstrj as the main file, then load bonds.dump as a topology overlay. The bond dump uses batom1/batom2/btype columns that OVITO recognises directly. For more detail please see the OVITO manual.

Verifying the output¶

Read back the final snapshot and check that crosslinks actually formed. The reset_mol_ids command at the end of the script reassigns molecule IDs based on bond connectivity — if the count dropped from 36, networks have formed.

In [13]:

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)")

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)")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[13], line 1
----> 1 final = mp.io.read_lammps_data(output_dir / "final.data", atom_style="full")
      2 n_atoms_final = final["atoms"].nrows
      3 n_mols_final = len(set(final["atoms"]["mol_id"]))
      4 print(f"final: {n_atoms_final} atoms, {n_mols_final} molecules (started with 36)")

NameError: name 'mp' is not defined

Troubleshooting¶

Template generation fails

Print the selected site and leaving-group atoms to see what the selectors found:

In [14]:

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')}")

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')}")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[14], line 1
----> 1 site = find_port(left, "$")
      2 carbon = [a for a in find_neighbors(left, site) if a.get("element") == "C"][0]
      3 print(f"site: {carbon.get('element')} name={carbon.get('name')}")
      4 for nb in find_neighbors(left, carbon):

NameError: name 'find_port' is not defined

Packing fails

Reduce the target density or monomer count. Packmol needs enough room to place molecules without overlap.

LAMMPS: "Atom type affected by reaction is too close to template edge"

Increase the radius parameter in BondReactReacter. A larger radius captures more of the local environment, pushing the template boundary further from type-changed atoms.

LAMMPS: reactions fire but molecule count stays constant

The post template is missing the new crosslinking bond. This happens when the typifier is not passed to reacter.run() — the new bond has no force field type and gets silently dropped during export. Always pass typifier=typifier.

LAMMPS rejects templates

Check that the .map file contains valid equivalence IDs:

In [15]:

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

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

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[15], line 1
----> 1 print((output_dir / "rxn1.map").read_text()[:500])

File ~/.asdf/installs/python/3.13.3/lib/python3.13/pathlib/_local.py:546, in Path.read_text(self, encoding, errors, newline)
    543 # Call io.text_encoding() here to ensure any warning is raised at an
    544 # appropriate stack level.
    545 encoding = io.text_encoding(encoding)
--> 546 return PathBase.read_text(self, encoding, errors, newline)

File ~/.asdf/installs/python/3.13.3/lib/python3.13/pathlib/_abc.py:632, in PathBase.read_text(self, encoding, errors, newline)
    628 def read_text(self, encoding=None, errors=None, newline=None):
    629     """
    630     Open the file in text mode, read it, and close the file.
    631     """
--> 632     with self.open(mode='r', encoding=encoding, errors=errors, newline=newline) as f:
    633         return f.read()

File ~/.asdf/installs/python/3.13.3/lib/python3.13/pathlib/_local.py:537, in Path.open(self, mode, buffering, encoding, errors, newline)
    535 if "b" not in mode:
    536     encoding = io.text_encoding(encoding)
--> 537 return io.open(self, mode, buffering, encoding, errors, newline)

FileNotFoundError: [Errno 2] No such file or directory: '04_output/rxn1.map'

See also: Topology-Driven Assembly, Polydisperse Systems.