Stepwise Polymer Building¶

"How do I connect monomers into a polymer — one observable step at a time?"

When you start working with polymers, the hardest part is often debuggability: did the monomers get placed correctly, did the reaction remove the right atoms, did the new bond form where you expected?

This guide demonstrates a minimal, auditable workflow: build a monomer with explicit reactive ports, place a second monomer in a controlled way, then connect them using a simple dehydration reaction. After each stage we export a LAMMPS file so you can inspect the structure in OVITO/VMD.

What We Will Do¶

1. Build an ethylene-oxide monomer from BigSMILES (ports define connection sites).
2. Create a second monomer and place it near the first (rotate + translate).
3. Connect the two fragments via a dehydration reaction (remove H₂O, form an ether bond).
4. Export checkpoints as LAMMPS data files for visualization.

Requirements & Outputs¶

• RDKit is required for 3D coordinate generation in this notebook (RDKitAdapter / Generate3D).
• Outputs are written to user-guide-output/02_polymer_stepwise/ (safe to delete).

If you only want to understand the chemistry + API flow and do not care about 3D coordinates, you can still read the notebook top-to-bottom; execution will fail at the RDKit step if RDKit is missing.

In [1]:

from pathlib import Path

import numpy as np

import molpy as mp
from molpy.core.atomistic import Atomistic
from molpy.adapter import RDKitAdapter
from molpy.compute import Generate3D
from molpy.io.data.lammps import LammpsDataWriter
from molpy.parser.smiles import bigsmilesir_to_monomer, parse_bigsmiles
from molpy.reacter import (
    Reacter,
    find_port_atom,
    form_single_bond,
    select_c_neighbor,
    select_hydroxyl_group,
    select_hydroxyl_h_only,
    select_port
    )
from molpy.typifier.atomistic import OplsAtomisticTypifier

from pathlib import Path import numpy as np import molpy as mp from molpy.core.atomistic import Atomistic from molpy.adapter import RDKitAdapter from molpy.compute import Generate3D from molpy.io.data.lammps import LammpsDataWriter from molpy.parser.smiles import bigsmilesir_to_monomer, parse_bigsmiles from molpy.reacter import ( Reacter, find_port_atom, form_single_bond, select_c_neighbor, select_hydroxyl_group, select_hydroxyl_h_only, select_port ) from molpy.typifier.atomistic import OplsAtomisticTypifier

Helper Functions¶

This notebook is easiest to follow if we keep two operations explicit and repeatable:

• Build a monomer from a BigSMILES string with ports ([<], [>]).
• Export a checkpoint to LAMMPS so you can inspect the exact geometry at that stage.

We deliberately keep the helpers small and “boring”: if something looks wrong later, you can debug it without digging through abstractions.

In [2]:

def build_monomer() -> Atomistic:
    """Build monomer with -OH end groups from BigSMILES."""
    bigsmiles = "{[][<]OCCOCCOCCO[>][]}"
    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()
    return monomer

def build_monomer() -> Atomistic: """Build monomer with -OH end groups from BigSMILES.""" bigsmiles = "{[][<]OCCOCCOCCO[>][]}" 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() return monomer

In [3]:

def export_frame_to_lammps(
    atomistic: Atomistic,
    output_path: Path,
    box_size: float = 20.0,
) -> None:
    """Export Atomistic to LAMMPS after typification."""
    frame = atomistic.to_frame()
    frame.metadata["box"] = mp.Box.cubic(length=box_size)

    n_atoms = frame["atoms"].nrows

    # Ensure ID field
    if "id" not in frame["atoms"]:
        ids = [atom.get("id", i) for i, atom in enumerate(atomistic.atoms, start=1)]
        frame["atoms"]["id"] = np.array(ids, dtype=int)

    # Ensure mol field
    frame["atoms"]["mol"] = np.array([1] * n_atoms, dtype=int)

    # Ensure charge field is available as 'q' (LAMMPS uses 'q')
    if "charge" in frame["atoms"]:
        frame["atoms"]["q"] = frame["atoms"]["charge"]

    writer = LammpsDataWriter(output_path, atom_style="full")
    writer.write(frame)

def export_frame_to_lammps( atomistic: Atomistic, output_path: Path, box_size: float = 20.0, ) -> None: """Export Atomistic to LAMMPS after typification.""" frame = atomistic.to_frame() frame.metadata["box"] = mp.Box.cubic(length=box_size) n_atoms = frame["atoms"].nrows # Ensure ID field if "id" not in frame["atoms"]: ids = [atom.get("id", i) for i, atom in enumerate(atomistic.atoms, start=1)] frame["atoms"]["id"] = np.array(ids, dtype=int) # Ensure mol field frame["atoms"]["mol"] = np.array([1] * n_atoms, dtype=int) # Ensure charge field is available as 'q' (LAMMPS uses 'q') if "charge" in frame["atoms"]: frame["atoms"]["q"] = frame["atoms"]["charge"] writer = LammpsDataWriter(output_path, atom_style="full") writer.write(frame)

Step 1: Load ForceField and Create Typifier¶

A Typifier is responsible for mapping atoms in a molecular structure to specific atom types defined in a force field. It examines the chemical environment of each atom—including its element and connectivity—to assign the appropriate parameters required for molecular dynamics simulations.

In [4]:

# Create output directory (safe to delete)
output_dir = Path("user-guide-output/02_polymer_stepwise")
output_dir.mkdir(parents=True, exist_ok=True)

# Load force field (built-in)
ff = mp.io.read_xml_forcefield("oplsaa.xml")
typifier = OplsAtomisticTypifier(ff, strict_typing=True)
print("Force field loaded.")
print(f"Outputs: {output_dir.resolve()}")

# Create output directory (safe to delete) output_dir = Path("user-guide-output/02_polymer_stepwise") output_dir.mkdir(parents=True, exist_ok=True) # Load force field (built-in) ff = mp.io.read_xml_forcefield("oplsaa.xml") typifier = OplsAtomisticTypifier(ff, strict_typing=True) print("Force field loaded.") print(f"Outputs: {output_dir.resolve()}")

2025-12-30 16:08:38,893 - 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.
Outputs: /opt/buildhome/repo/docs/user-guide/user-guide-output/02_polymer_stepwise

Step 2: Build and Export First Monomer¶

Process:

1. Build monomer from BigSMILES
2. Generate topology (bonds, angles, dihedrals)
3. Assign force field types
4. Merge into polymer container
5. Export to LAMMPS

Why generate topology? This topology is required by classical molecular dynamics to define bonded interaction terms.

Why merge into polymer? This creates a container that will hold the growing polymer chain.

In [5]:

monomer = build_monomer()
monomer.get_topo(gen_angle=True, gen_dihe=True)
typifier.typify(monomer)

polymer = Atomistic()
polymer.merge(monomer)

print("Step 4 completed: first monomer built and typed")
print(f"  Monomer: {len(monomer.atoms)} atoms")
print(f"  Polymer: {len(polymer.atoms)} atoms")

export_frame_to_lammps(polymer, output_dir / "step1_monomer.data", box_size=20.0)
print(f"  Exported: {output_dir / 'step1_monomer.data'}")

left_monomer = monomer

monomer = build_monomer() monomer.get_topo(gen_angle=True, gen_dihe=True) typifier.typify(monomer) polymer = Atomistic() polymer.merge(monomer) print("Step 4 completed: first monomer built and typed") print(f" Monomer: {len(monomer.atoms)} atoms") print(f" Polymer: {len(polymer.atoms)} atoms") export_frame_to_lammps(polymer, output_dir / "step1_monomer.data", box_size=20.0) print(f" Exported: {output_dir / 'step1_monomer.data'}") left_monomer = monomer

Step 4 completed: first monomer built and typed
  Monomer: 24 atoms
  Polymer: 24 atoms
  Exported: user-guide-output/02_polymer_stepwise/step1_monomer.data

Step 3: Add Second Monomer with Positioning¶

We place the second monomer using only two operations:

1. Rotate (to face the ports)
2. Translate (move it to the target location)

For this example, we use a fixed separation distance to avoid overlaps.

In [6]:

monomer2 = build_monomer()

# Assign new IDs (must be unique after merge)
max_id = max((atom.get("id", 0) for atom in polymer.atoms), default=0)
for atom in monomer2.atoms:
    max_id += 1
    atom["id"] = max_id

# Ports: head-to-head uses ">" on both fragments
port1_atom = find_port_atom(polymer, ">")
port2_atom = find_port_atom(monomer2, ">")

port1_pos = np.array([port1_atom["x"], port1_atom["y"], port1_atom["z"]], dtype=float)
port2_pos = np.array([port2_atom["x"], port2_atom["y"], port2_atom["z"]], dtype=float)

polymer_center = np.mean(np.array([[a["x"], a["y"], a["z"]] for a in polymer.atoms], dtype=float), axis=0)
monomer2_center = np.mean(np.array([[a["x"], a["y"], a["z"]] for a in monomer2.atoms], dtype=float), axis=0)

direction = (port1_pos - polymer_center) / np.linalg.norm(port1_pos - polymer_center)
port2_direction = (port2_pos - monomer2_center) / np.linalg.norm(port2_pos - monomer2_center)

# Rotate
axis = np.cross(port2_direction, -direction)
monomer2.rotate(axis=axis.tolist(), angle=np.pi, about=monomer2_center.tolist())

# Translate (hardcoded separation)
separation = 10.0  # Å
port2_pos = np.array([port2_atom["x"], port2_atom["y"], port2_atom["z"]], dtype=float)
target_pos = port1_pos + direction * separation
monomer2.move(delta=(target_pos - port2_pos).tolist())

# Typify and merge
monomer2.get_topo(gen_angle=True, gen_dihe=True)
typifier.typify(monomer2)

polymer.merge(monomer2)
polymer.get_topo(gen_angle=True, gen_dihe=True)
typifier.typify(polymer)

print("Step 5 completed: second monomer rotated + translated, merged, and typed")
export_frame_to_lammps(polymer, output_dir / "step2_two_monomers.data", box_size=25.0)
print(f"  Exported: {output_dir / 'step2_two_monomers.data'}")

left_monomer_copy = left_monomer.copy()
right_monomer_copy = monomer2.copy()

monomer2 = build_monomer() # Assign new IDs (must be unique after merge) max_id = max((atom.get("id", 0) for atom in polymer.atoms), default=0) for atom in monomer2.atoms: max_id += 1 atom["id"] = max_id # Ports: head-to-head uses ">" on both fragments port1_atom = find_port_atom(polymer, ">") port2_atom = find_port_atom(monomer2, ">") port1_pos = np.array([port1_atom["x"], port1_atom["y"], port1_atom["z"]], dtype=float) port2_pos = np.array([port2_atom["x"], port2_atom["y"], port2_atom["z"]], dtype=float) polymer_center = np.mean(np.array([[a["x"], a["y"], a["z"]] for a in polymer.atoms], dtype=float), axis=0) monomer2_center = np.mean(np.array([[a["x"], a["y"], a["z"]] for a in monomer2.atoms], dtype=float), axis=0) direction = (port1_pos - polymer_center) / np.linalg.norm(port1_pos - polymer_center) port2_direction = (port2_pos - monomer2_center) / np.linalg.norm(port2_pos - monomer2_center) # Rotate axis = np.cross(port2_direction, -direction) monomer2.rotate(axis=axis.tolist(), angle=np.pi, about=monomer2_center.tolist()) # Translate (hardcoded separation) separation = 10.0 # Å port2_pos = np.array([port2_atom["x"], port2_atom["y"], port2_atom["z"]], dtype=float) target_pos = port1_pos + direction * separation monomer2.move(delta=(target_pos - port2_pos).tolist()) # Typify and merge monomer2.get_topo(gen_angle=True, gen_dihe=True) typifier.typify(monomer2) polymer.merge(monomer2) polymer.get_topo(gen_angle=True, gen_dihe=True) typifier.typify(polymer) print("Step 5 completed: second monomer rotated + translated, merged, and typed") export_frame_to_lammps(polymer, output_dir / "step2_two_monomers.data", box_size=25.0) print(f" Exported: {output_dir / 'step2_two_monomers.data'}") left_monomer_copy = left_monomer.copy() right_monomer_copy = monomer2.copy()

Step 5 completed: second monomer rotated + translated, merged, and typed
  Exported: user-guide-output/02_polymer_stepwise/step2_two_monomers.data

Step 4: Connect Monomers via Dehydration Reaction¶

Reaction Mechanism¶

Dehydration reaction: R-OH + HO-R' → R-O-R' + H₂O

Reacter Configuration:

• anchor_selector_left: Selects C neighbor of -OH (where new bond forms)
• anchor_selector_right: Selects O atom itself
• leaving_selector_left: Selects entire -OH group
• leaving_selector_right: Selects only H from -OH
• bond_former: Creates C-O single bond

Why this design?

• Left side: C-OH → C (remove OH, keep C)
• Right side: HO-R → O-R (remove H, keep O)
• Result: C-O bond formation, H₂O removed

Port atoms: We use find_port_atom() to locate the reactive sites marked in BigSMILES.

In [7]:

# Define dehydration reaction
dehydration = Reacter(
    name="dehydration_ether",
    anchor_selector_left=select_c_neighbor,  # C next to -OH
    anchor_selector_right=select_port,  # O atom itself
    leaving_selector_left=select_hydroxyl_group,  # remove -OH
    leaving_selector_right=select_hydroxyl_h_only,  # remove H only
    bond_former=form_single_bond,  # create C-O bond
)

# Run reaction
result = dehydration.run(
    left=left_monomer_copy,
    right=right_monomer_copy,
    port_atom_L=find_port_atom(left_monomer_copy, ">"),
    port_atom_R=find_port_atom(right_monomer_copy, ">"),
    compute_topology=True,
)

product = result.product_info.product

# Typify product
product.get_topo(gen_angle=True, gen_dihe=True)
typifier.typify(product)

print("Step 6 completed: monomers connected via dehydration")
print(f"  Reacted polymer: {len(product.atoms)} atoms")
print(f"  Removed atoms: {len(result.topology_changes.removed_atoms)} atoms (water: O+H+H)")
print(f"  New bonds: {len(result.topology_changes.new_bonds)} bonds")

export_frame_to_lammps(product, output_dir / "step3_connected.data", box_size=25.0)
print(f"  Exported: {output_dir / 'step3_connected.data'}")

# Define dehydration reaction dehydration = Reacter( name="dehydration_ether", anchor_selector_left=select_c_neighbor, # C next to -OH anchor_selector_right=select_port, # O atom itself leaving_selector_left=select_hydroxyl_group, # remove -OH leaving_selector_right=select_hydroxyl_h_only, # remove H only bond_former=form_single_bond, # create C-O bond ) # Run reaction result = dehydration.run( left=left_monomer_copy, right=right_monomer_copy, port_atom_L=find_port_atom(left_monomer_copy, ">"), port_atom_R=find_port_atom(right_monomer_copy, ">"), compute_topology=True, ) product = result.product_info.product # Typify product product.get_topo(gen_angle=True, gen_dihe=True) typifier.typify(product) print("Step 6 completed: monomers connected via dehydration") print(f" Reacted polymer: {len(product.atoms)} atoms") print(f" Removed atoms: {len(result.topology_changes.removed_atoms)} atoms (water: O+H+H)") print(f" New bonds: {len(result.topology_changes.new_bonds)} bonds") export_frame_to_lammps(product, output_dir / "step3_connected.data", box_size=25.0) print(f" Exported: {output_dir / 'step3_connected.data'}")

Step 6 completed: monomers connected via dehydration
  Reacted polymer: 45 atoms
  Removed atoms: 3 atoms (water: O+H+H)
  New bonds: 1 bonds
  Exported: user-guide-output/02_polymer_stepwise/step3_connected.data

Summary¶

We demonstrated a step-by-step polymer construction workflow:

1. Built a monomer from BigSMILES (ports define reactive sites)
2. Positioned a second monomer with explicit translation + rotation
3. Connected the fragments using a Reacter dehydration mechanism

Key Takeaways

• BigSMILES makes reactive sites explicit and auditable.
• Geometry placement is a separate concern from chemistry; keep it deterministic.
• After any edit (merge/reaction), regenerate topology and re-run typing.

Visualization Open the exported .data files in OVITO/VMD to compare steps 1–3.