Crosslinked Networks¶
After reading this page you will be able to generate fix bond/react templates, pack a crosslinkable monomer mixture, and export all files LAMMPS needs for reactive MD.
!!! note "Prerequisites"
This guide requires RDKit, Packmol, and the oplsaa.xml force field. Familiarity with Stepwise Polymer Construction is assumed.
A crosslinking simulation needs four deliverables¶
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.
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)
2026-04-10 07:26:48,934 - 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.
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.
eo2 = mp.parser.parse_monomer("{[][$]OCCOCCOCCO[$][]}")
eo2 = mp.tool.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)
eo3 = mp.parser.parse_monomer("{[]C(COCCO[$])(COCCO[$])COCCO[$][]}")
eo3 = mp.tool.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)
/opt/buildhome/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/tool/rdkit.py:110: 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)
Verify that the port markers survived 3D generation:
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}")
EO2: atoms=24, ports=['$', '$'] EO3: atoms=38, ports=['$', '$', '$']
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.
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")
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.
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
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.
print(f"pre: {len(template.pre.atoms)} atoms")
print(f"post: {len(template.post.atoms)} atoms")
pre: 23 atoms post: 23 atoms
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.
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)
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")
packed: 990 atoms in 21.1 Å box
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.
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},
)
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.
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.
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:]}"
--------------------------------------------------------------------------- FileNotFoundError Traceback (most recent call last) Cell In[12], line 7 3 script = mp.Script.from_text("run", lammps_script, language="other") 4 script.tags.add("input") 5 6 engine = LAMMPSEngine(executable="lmp_serial", check_executable=False) ----> 7 proc = engine.run( 8 script, 9 workdir=output_dir, 10 capture_output=True, File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/engine/base.py:274, in Engine.run(self, scripts, workdir, capture_output, check, timeout, **kwargs) 270 script.save(script_path) 272 self.input_script = self._find_input_script() --> 274 return self._execute( 275 run_dir, 276 capture_output=capture_output, 277 check=check, 278 timeout=timeout, 279 **kwargs, 280 ) File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/engine/lammps.py:115, in LAMMPSEngine._execute(self, run_dir, capture_output, check, timeout, **kwargs) 110 input_file = self.input_script.path.name 111 command = self._build_full_command( 112 ["-in", input_file, "-log", "log.lammps", "-screen", "none"] 113 ) --> 115 return subprocess.run( 116 command, 117 cwd=run_dir, 118 capture_output=capture_output, 119 text=True, 120 check=check, 121 timeout=timeout, 122 env=self._merged_env(), 123 ) File ~/.asdf/installs/python/3.13.3/lib/python3.13/subprocess.py:554, in run(input, capture_output, timeout, check, *popenargs, **kwargs) 551 kwargs['stdout'] = PIPE 552 kwargs['stderr'] = PIPE --> 554 with Popen(*popenargs, **kwargs) as process: 555 try: 556 stdout, stderr = process.communicate(input, timeout=timeout) File ~/.asdf/installs/python/3.13.3/lib/python3.13/subprocess.py:1039, in Popen.__init__(self, args, bufsize, executable, stdin, stdout, stderr, preexec_fn, close_fds, shell, cwd, env, universal_newlines, startupinfo, creationflags, restore_signals, start_new_session, pass_fds, user, group, extra_groups, encoding, errors, text, umask, pipesize, process_group) 1035 if self.text_mode: 1036 self.stderr = io.TextIOWrapper(self.stderr, 1037 encoding=encoding, errors=errors) -> 1039 self._execute_child(args, executable, preexec_fn, close_fds, 1040 pass_fds, cwd, env, 1041 startupinfo, creationflags, shell, 1042 p2cread, p2cwrite, 1043 c2pread, c2pwrite, 1044 errread, errwrite, 1045 restore_signals, 1046 gid, gids, uid, umask, 1047 start_new_session, process_group) 1048 except: 1049 # Cleanup if the child failed starting. 1050 for f in filter(None, (self.stdin, self.stdout, self.stderr)): File ~/.asdf/installs/python/3.13.3/lib/python3.13/subprocess.py:1969, in Popen._execute_child(self, args, executable, preexec_fn, close_fds, pass_fds, cwd, env, startupinfo, creationflags, shell, p2cread, p2cwrite, c2pread, c2pwrite, errread, errwrite, restore_signals, gid, gids, uid, umask, start_new_session, process_group) 1967 err_msg = os.strerror(errno_num) 1968 if err_filename is not None: -> 1969 raise child_exception_type(errno_num, err_msg, err_filename) 1970 else: 1971 raise child_exception_type(errno_num, err_msg) FileNotFoundError: [Errno 2] No such file or directory: 'lmp_serial'
!!! 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.
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)")
--------------------------------------------------------------------------- FileNotFoundError 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)") File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/io/readers.py:44, in read_lammps_data(file, atom_style, frame) 41 from .data.lammps import LammpsDataReader 43 reader = LammpsDataReader(Path(file), atom_style) ---> 44 return reader.read(frame=frame) File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/io/data/lammps.py:51, in LammpsDataReader.read(self, frame) 48 frame = frame or Frame() 50 # Read and parse the file ---> 51 lines = self._read_lines() 52 sections = self._extract_sections(lines) 54 # Parse header and set up box File ~/.asdf/installs/python/3.13.3/lib/python3.13/site-packages/molpy/io/data/lammps.py:126, in LammpsDataReader._read_lines(self) 124 def _read_lines(self) -> list[str]: 125 """Read file and return non-empty, non-comment lines.""" --> 126 with open(self._path) as f: 127 return [ 128 line.strip() 129 for line in f 130 if line.strip() and not line.strip().startswith("#") 131 ] FileNotFoundError: [Errno 2] No such file or directory: '04_output/final.data'
Troubleshooting¶
Template generation fails
Print the selected site and leaving-group atoms to see what the selectors found:
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: C name=None neighbor: O name=None neighbor: C name=None neighbor: H name=None neighbor: H name=None
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:
print((output_dir / "rxn1.map").read_text()[:500])
# auto-generated map file for fix bond/react 23 equivalences 2 edgeIDs 3 deleteIDs InitiatorIDs 18 2 EdgeIDs 6 16 DeleteIDs 1 7 15 Equivalences 1 1 2 2 3 3 4 4 5 5 6 6 7 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 21 21 22 22 23 23
See also: Topology-Driven Assembly, Polydisperse Systems.