Builder¶

System assembly: select every reaction, execute the graph edits as one batch, retype the region each edit disturbed, then optionally finalize topology. Growing a chain and crosslinking a melt are the same algorithm with a different pairing rule, so there is one kernel and one variation point.

Quick reference¶

Canonical example¶

A repeat unit is an ordinary capped molecule with a few of its atoms named. There is no port system and no direction: the reaction SMARTS is the only place the chemistry lives, and %a / %b bind it to the atoms you marked.

import molpy as mp
from molpy.builder.assembly import (
    MonomerLibrary,
    PolymerBuilder,
    ResiduePlacer,
    SiteMap,
)
from molpy.conformer import Conformer
from molpy.core import fields

eo, _ = Conformer(add_hydrogens=True, seed=42).generate(
    mp.io.read_smiles("OCCO")
)
SiteMap(eo).label_elements("O", "a", "b")

ether = mp.Reaction("[O;%a:1][H].[C:2][O;%b][H]>>[O:1][C:2]")
builder = PolymerBuilder(MonomerLibrary({"EO": eo}), ether, placer=ResiduePlacer())
chain = builder.build_linear("EO", 5)

assert chain.__class__.__name__ == "Atomistic"
assert len({int(a[fields.RES_ID]) for a in chain.atoms}) == 5

Each repeat unit is a residue, and that identity survives into a PDB or a prmtop — it is output, not a build-time marker to scrub afterwards.

React first; finalize when needed¶

typifier= never receives the growing polymer. The builder executes all selected reactions in one batch against the intact world, then cuts an AffectedRegion around each edit and types that. The order matters: a region cut before the edit would have to model the edit itself — including every neighbouring reaction close enough to change a type, and the unmapped atoms each one consumes — while a region cut afterwards inherits a graph on which those edits already succeeded.

Structurally identical junctions share one typing pass via RetypeCache, and only the region's interior is written back. Bond/angle/dihedral annotations are never copied out of a region.

The default finalization generates complete angle/dihedral topology once. For a large system that will be written to an MD format later, select Finalization.ATOMS while building and apply StructureFinalizer(Finalization.TOPOLOGY) at export time. This keeps the public examples focused on polymer architectures while still documenting the deferred topology option.

Use Finalization.BONDED together with bonded=ForceFieldParams(forcefield) when those generated terms also need type and parameter annotations.

A full walk of architectures from one ethylene-glycol template — linear, ring, star, comb, gels, dual network — is the user-guide section Polymer Topologies (paired with examples/topology/).

Nanostructure topology¶

Nanostructure builders keep their lattice planning private and expose two products: the molecular graph, and the simulation cell it was laid out in.

from molpy.builder import CarbonTubeBuilder

tube_builder = CarbonTubeBuilder(6, 6, cells=2, periodic=True)
tube = tube_builder.build()
assert len(tube.atoms) == 48
assert len(tube.bonds) == 72
assert tube_builder.cell().pbc.tolist() == [False, False, True]

The same builder accepts zigzag (n, 0), armchair (n, n), and general chiral (n, m) tubes. See Nanostructures for open ends, length selection, and deferred topology.

Crosslinking is the same machine¶

Strip the library and the notation away and you have the kernel itself, which is all crosslinking needs: a graph you already have, plus a rule for which sites pair up.

from molpy.builder.assembly import GraphAssembler, RandomSelector

melt = mp.Atomistic()
for i in range(4):
    melt.def_atom(element="N", x=float(i), y=0.0, z=0.0)
    melt.def_atom(element="O", x=float(i), y=1.0, z=0.0)

gel = GraphAssembler(mp.Reaction("[N:1].[O:2]>>[N:1][O:2]")).assemble(
    melt, RandomSelector(conversion=1.0, seed=1, cutoff=2.0)
)
assert len(list(gel.bonds)) == 4

No placer is passed here: the melt's coordinates are already meaningful and must not be disturbed. That is a decision about your input, not about which class you reached for, which is why it is an argument.

Polydisperse systems¶

Sample a chain plan, then loop build:

import numpy as np
from molpy.builder.polymer import (
    PolydisperseChainGenerator,
    SchulzZimmPolydisperse,
    SystemPlanner,
    WeightedSequenceGenerator,
)

planner = SystemPlanner(
    PolydisperseChainGenerator(
        WeightedSequenceGenerator({"EO": 1.0}),
        {"EO": 44.05},
        distribution=SchulzZimmPolydisperse(1500, 3000),
    ),
    target_total_mass=5e3,
)
plan = planner.plan_system(np.random.default_rng(42))
chains = [builder.build_linear("EO", len(c.monomers)) for c in plan.chains[:2]]
assert len(chains) == 2

Your own pairing rule¶

You never subclass the assembler. You write a Selector, which answers the one question that varies. The matching has already happened — the kernel does it once, in linear time — so a selector never scans the system, it only decides.

from molpy.builder.assembly import Selector


class FirstPairSelector(Selector):
    """React exactly one pairing: the first site of each reactant."""

    def select(self, context):
        a_sites = context.occurrences[context.comp_a]
        b_sites = context.occurrences[context.comp_b]
        yield {**a_sites[0], **b_sites[0]}


one = GraphAssembler(mp.Reaction("[N:1].[O:2]>>[N:1][O:2]")).assemble(
    melt, FirstPairSelector()
)
assert len(list(one.bonds)) == 1

Related¶

• Guide: Assembly