Atomistic and Topology¶

How do you add a leaving group, break a bond, or ask “is this carbon aromatic” before any coordinates go near an engine?

A coordinate list cannot answer. You need a connectivity model: which atoms are bonded, and what follows from that. In MolPy that model is an editable molecular graph.

Atomistic is a molecular graph you edit. Atoms are nodes, bonds are edges. Angles, dihedrals, and rings are derived from connectivity when you ask — not hand-maintained tables that go stale.

What it is not: a simulation snapshot, a force-field file, or a packed box. Those live on Frame and friends after the chemistry is settled (Block and Frame).

Why a graph?¶

Simulation tools face two kinds of question.

Geometric questions — “how far apart are these two atoms?” — need coordinates. Chemical questions — “which atoms share a covalent bond?”, “what happens when I break this bond?” — need connectivity.

A plain list of positions answers only the first kind. A molecular graph answers both identities and relations: you know which carbon is bonded to which oxygen, not merely that row 3 sits near row 5. That is why building, editing, and reacting molecules happen here, before the system is handed to an engine.

Building a molecule¶

Start empty, then add atoms and bonds. Properties are free-form keywords — element, coordinates, charge, or anything else your workflow needs.

import molpy as mp

mol = mp.Atomistic(name="ethanol")
c1 = mol.def_atom(element="C", name="C1", x=0.0, y=0.0, z=0.0)
c2 = mol.def_atom(element="C", name="C2", x=1.54, y=0.0, z=0.0)
o = mol.def_atom(element="O", name="O1", x=2.0, y=1.4, z=0.0)
h_o = mol.def_atom(element="H", name="HO", x=2.9, y=1.4, z=0.0)

mol.def_bond(c1, c2, order=1)
mol.def_bond(c2, o, order=1)
mol.def_bond(o, h_o, order=1)
print(f"{len(mol.atoms)} atoms, {len(mol.bonds)} bonds")
# -> 4 atoms, 3 bonds

Those atoms and bonds are live objects inside the graph — not detached copies. Change a property and the molecule changes immediately.

print(c1["element"])  # -> C
c1["charge"] = -0.18
print(c1["charge"])  # -> -0.18

bond = mol.bonds[0]
print(bond.itom["name"], bond.jtom["name"], bond.get("order"))
# -> C1 C2 1

Connectivity lives in the molecule, not in the atom¶

An atom does not store its own neighbour list. The container owns connectivity, so graph operations stay explicit and consistent. topo_distances walks the bond graph from one atom handle and returns (handle, hops) pairs; the neighbours are the atoms one bond away.

neighbors = [h for h, d in mol.topo_distances(c2.handle, max_hops=1) if d == 1]
print(sorted(mol.get(h, "name") for h in neighbors))
# -> ['C1', 'O1']

Removing an atom removes incident bonds with it — you never keep a dangling edge.

print(f"Before: {len(mol.atoms)} atoms, {len(mol.bonds)} bonds")
# -> Before: 4 atoms, 3 bonds
mol.del_atom(h_o)
print(f"After:  {len(mol.atoms)} atoms, {len(mol.bonds)} bonds")
# -> After:  3 atoms, 2 bonds

Topology is derived from bonds¶

Engines need angles and dihedrals as well as bonds. Maintaining those lists by hand is brittle: every bond edit would require a matching update.

MolPy treats bonded topology as derived from the bond graph. generate_topology reads the current bonds and writes perceived angles and dihedrals into the same Atomistic (in place; it returns the counts added, not the molecule). If the bond graph changes later, call generate_topology again.

propane = mp.Atomistic(name="propane")
ca = propane.def_atom(element="C", name="C1", x=0.0, y=0.0, z=0.0)
cb = propane.def_atom(element="C", name="C2", x=1.54, y=0.0, z=0.0)
cc = propane.def_atom(element="C", name="C3", x=3.08, y=0.0, z=0.0)
propane.def_bond(ca, cb)
propane.def_bond(cb, cc)

print(len(propane.angles), len(propane.dihedrals))
# -> 0 0
propane.generate_topology(gen_angle=True, gen_dihedral=True)
print(len(propane.angles), len(propane.dihedrals))
# -> 1 0

for angle in propane.angles:
    print(" — ".join(a["name"] for a in angle.endpoints))
# -> C1 — C2 — C3

Graph-distance queries (neighbours within \(n\) bonds, BFS distances) run on the same structure once connectivity is defined:

within_one = propane.topo_distances(cb.handle, max_hops=1)
print(sorted(propane.get(h, "name") for h, _ in within_one))
# -> ['C1', 'C2', 'C3']
dists = propane.topo_distances(ca.handle)
print({propane.get(h, "name"): d for h, d in dists})
# -> {'C1': 0, 'C2': 1, 'C3': 2}

Composition and copies¶

Independent clones use copy(). merge adds another graph into this one in place (it returns the old-to-new handle map); many copies are repeated merges of translated copies.

water = mp.Atomistic(name="water")
ow = water.def_atom(element="O", x=0.0, y=0.0, z=0.0)
h1 = water.def_atom(element="H", x=0.957, y=0.0, z=0.0)
h2 = water.def_atom(element="H", x=-0.239, y=0.927, z=0.0)
water.def_bond(ow, h1)
water.def_bond(ow, h2)

two = water.copy()
two.merge(water.copy().translate([5.0, 0.0, 0.0]))
print(len(two.atoms), len(two.bonds))
# -> 6 4

row = mp.Atomistic(name="row")
for i in range(4):
    row.merge(water.copy().translate([i * 4.0, 0.0, 0.0]))
print(len(row.atoms))
# -> 12

Bulk numeric access

When you need every \(x\) coordinate as an array — not each atom as a Python object — use mol.atoms["x"], mol.atoms["x", "y", "z"], or mol.column("x"). Prefer list(mol.atoms) only when you need identity-stable objects for editing or graph algorithms.

When to stay here, when to leave¶

Stay on Atomistic while the chemistry is still under discussion: add atoms, define bonds, inspect connectivity, run reactions.

When the chemistry is stable and the next job is export, analysis, or simulation, move to arrays: Block and Frame.

Check yourself¶

Why can two atoms with identical element and coordinates still be different atoms in MolPy?
After del_atom on a terminal hydrogen, how many bonds should remain on the heavy-atom skeleton of the ethanol example above?
If you add a bond after generate_topology, are the angle lists automatically up to date? What do you call?