Molecules and geometry¶

Matmalib represents a geometry as a State, containing an ordered list of Atom objects. A Fragment groups existing atoms, and a Space collects states. Coordinates are ordinary NumPy arrays, so the same objects work with the conformer and carbohydrate tools.

This tutorial builds a molecule, measures internal coordinates, edits an acyclic chain, and separates disconnected fragments. See core reference notes for the complete method inventory and current limitations.

Build a state from coordinates¶

Always provide three coordinates to Atom, use floating-point coordinates when editing, and pass index= explicitly. The third positional argument is valence, which is currently unused. Use angstroms for molecular coordinates so the default bond thresholds have their intended meaning.

import numpy as np
from matmalib import Atom, State

water = State(
    id="water",
    atoms=[
        Atom(coord=[0.0, 0.0, 0.0], type="O", index=0),
        Atom(coord=[0.9572, 0.0, 0.0], type="H", index=1),
        Atom(coord=[-0.2390, 0.9270, 0.0], type="H", index=2),
    ],
    NAtoms=3,
)
for atom, mass in zip(water.atoms, [15.999, 1.008, 1.008]):
    atom.mass = mass

print("Coordinates (angstrom):\n", water.xyz)
print("Atomic masses (u):", water.Mass)
print("Total mass (u):", water.mass)
print("O-H distance:", round(water.measure_internal([0, 1]), 4))
angle, axis = water.measure_internal([1, 0, 2])
print("H-O-H angle (degrees):", round(angle, 2))

The O–H bond is 0.9572 Å; the angle is approximately 104.46°. xyz has shape (n_atoms, 3), Mass has shape (n_atoms,), and mass sums the supplied atomic masses. Matmalib does not populate atomic masses automatically. The example uses conventional average masses; use isotope masses if your calculation requires them.

Try changing: move the third atom, run again, and compare its angle. Masses do not affect these geometric measurements.

State(**kwargs) stores the attributes you supply; it does not validate or complete a molecular model. Supply atoms, NAtoms, and id when using connectivity, visualization, or writers. xyz and Mass assemble new arrays each time: edit state.atoms[i].coord, rather than expecting state.xyz[i] = ... to persist. The atom's x, y, and z attributes are snapshots created at initialization and do not track later coordinate edits; read coord instead.

Look up element symbols¶

State.element_symbol(8) returns "O"; StateBuilder.element_symbol("8") uses the same mapping. In the other direction, InputBuilder.element_number("O") returns 8. The underlying ELEMENT_SYMBOLS and ELEMENT_NUMBERS dictionaries cover elements 1–118, with string atomic-number keys and case-sensitive element symbols respectively. They provide identifiers, not atomic masses or valence information. The element utilities example lets you run all three wrappers and inspect the dictionaries.

Measure distances, angles, and dihedrals¶

Atom selections use zero-based positions in state.atoms. The older distance/angle editing and connectivity routines also require atom.id to equal that position. Dihedral measurements and edits use positions independently of Atom.id and accept Python or NumPy integers; other legacy geometry methods require Python int values. Do not mix integer positions and Atom objects in a selection.

import numpy as np
from matmalib import Atom, State, StateEditor

points = [[0.0, 1.0, 0.0], [0.0, 0.0, 0.0],
          [1.4, 0.0, 0.0], [1.4, 0.0, 1.0]]
chain = State(id="four_atoms", NAtoms=4, atoms=[
    Atom(coord=p, type="C", index=i) for i, p in enumerate(points)
])

print("Distance:", StateEditor.get_distance([1, 2], state=chain))
angle, axis = StateEditor.get_angle([0, 1, 2], state=chain)
print("Angle and axis:", angle, axis)
print("Dihedral:", StateEditor.get_dihedral([0, 1, 2, 3], state=chain))
print("Atom-object selection:", StateEditor.measure(chain.atoms[:2]))

for atom in chain.atoms:
    atom.coord += np.array([10.0, -4.0, 7.0])
print("Dihedral after translation:", chain.measure_internal([0, 1, 2, 3]))

This idealized geometry gives a 1.4 Å central distance, a 90° angle, and a +90° dihedral, which remains +90° after translation. The angle is centered on the middle atom. A dihedral describes the planes formed by atoms (0, 1, 2) and (1, 2, 3). Its positive direction follows a right-handed rotation about the central bond from atom 1 toward atom 2; see the signed torsion convention for the formula.

Selection	`State.measure_internal` / `StateEditor.measure` return
Two atoms	Distance as a float, in input coordinate units
Three atoms	`(angle_degrees, axis_vector)`; the axis is an unnormalized cross product
Four atoms	Signed dihedral in degrees

get_dihedral(..., ax=True) returns (dihedral_degrees, axis_vector). This compatibility axis is the unnormalized cross product of the two plane normals; it can be zero for a planar geometry, even though its dihedral is defined. get_distance, get_angle, and get_dihedral accept either all Atom objects, or all integer positions with a state. get_dihedral requires exactly four atoms with finite three-component coordinates and rejects negative/out-of-range positions, coincident central atoms, or collinear triples with ValueError. The older distance and angle length checks only print warnings. measure returns None for a selection length other than two, three, or four.

Try changing: change the last coordinate in points from 1.0 to -1.0. This reflection changes the dihedral to −90°, and translating the whole geometry still leaves it unchanged. For ring analysis use the dedicated puckering functions.

Compatibility: StateEditor.get_dihedral, StateEditor.measure, and State.measure_internal now use the same signed atan2 convention as the glycan torsion routines. This corrects the reversed legacy sign and its translation-dependent error. Existing glycan and ring-analysis results are unchanged. Recompute previously saved legacy dihedrals from coordinates when possible; simply negating them cannot repair values affected by the old translation bug.

Infer bonds and separate molecules¶

connectivity_matrix() sets three attributes: conn_mat, a symmetric zero/one matrix; graph, a NetworkX graph; and Nmols, the number of connected components. It returns None. A heavy-atom pair is connected when its distance is strictly less than distXX (default 1.65 Å), and any pair involving hydrogen uses distXH (default 1.15 Å).

from matmalib import Atom, State, StateEditor

points = [[0.0, 0.0, 0.0], [1.4, 0.0, 0.0],
          [5.0, 0.0, 0.0], [6.4, 0.0, 0.0]]
state = State(id="two_fragments", NAtoms=4, atoms=[
    Atom(coord=p, type="C", index=i) for i, p in enumerate(points)
])
state.connectivity_matrix()
state.system_split()

print("Bond matrix:\n", state.conn_mat)
print("Connected components:", state.Nmols)
print("Neighbors of atom 0:", StateEditor.adjacent_atoms(state.conn_mat, 0))
for fragment in state.mols:
    print("Fragment atom IDs:", [atom.id for atom in fragment.atoms])

You should find two fragments: atoms [0, 1] and [2, 3]. system_split() reads the existing state.graph and stores Fragment instances in state.mols. Fragments share Atom objects with their parent state, so moving a fragment atom also moves the corresponding state atom. Their xyz and Mass properties assemble the member coordinates and masses.

Try changing: set distXX=4.0 and see how the components merge. This illustrates why a distance cutoff is a bonding heuristic, not a bond-order or valence model. There are no periodic boundary conditions, minimum-image distances, or element-specific covalent radii in this method.

The legacy hydrogen ambiguity branch currently calls a nonexistent list method if a hydrogen is close to more than one possible neighbor. For problematic geometries, supply an explicit connectivity matrix and graph instead of relying on that branch. Carbohydrate auto-detection has its own connectivity implementation, described in glycans.

Edit an acyclic chain¶

Geometry edits move the connected component on the second side of a selected bond. Create connectivity first, work on a copy when preserving the input matters, and remeasure after an edit. These are geometric transformations, not energy minimization.

import copy
from matmalib import Atom, State, StateEditor

points = [[0.0, 1.0, 0.0], [0.0, 0.0, 0.0],
          [1.4, 0.0, 0.0], [1.4, 0.0, 1.0]]
original = State(id="chain", NAtoms=4, atoms=[
    Atom(coord=p, type="C", index=i) for i, p in enumerate(points)
])
original.connectivity_matrix()
edited = copy.deepcopy(original)
editor = StateEditor(edited)

carried = editor.determine_carried_atoms(edited, [1, 2])
print("Atoms moved with atom 2:", [atom.id for atom in carried])
editor.set_distance([1, 2], new_distance=1.5)
editor.set_angle([0, 1, 2], new_ang=110.0)
editor.set_dihedral([0, 1, 2, 3], new_dih=60.0)

print("Edited distance:", round(edited.measure_internal([1, 2]), 3))
print("Edited angle:", round(edited.measure_internal([0, 1, 2])[0], 3))
print("Edited dihedral:", round(edited.measure_internal([0, 1, 2, 3]), 3))
print("Original distance:", original.measure_internal([1, 2]))

The moved atoms are [2, 3]; the result measures 1.5 Å, 110°, and 60°, while the original retains its 1.4 Å central bond.

Try changing: choose a target dihedral of −60°, 0°, or 180°, or a distance of 1.6 Å. Re-run the full cell to start each experiment from the same geometry. With editor.set_dihedral(..., incr=True), new_dih is an increment in degrees rather than an absolute target: adding +30° to +60° gives +90°. Targets and increments wrap modulo 360°, so 780° is equivalent to 60°. Exact planar targets of 0° and ±180° are supported; +180° and −180° describe the same conformation.

Operation	Parameters and effect
`state.set_internal(atom_list, val)`	Dispatches by selection length to a distance, angle, or dihedral edit; mutates atoms and returns `None`
`editor.set_distance([a, b], new_distance)`	Translates the component on `b`'s side of bond `a–b`; distance is in input units
`editor.set_angle([a, b, c], new_ang)`	Rotates `c`'s component around `b`; target is in degrees
`editor.set_dihedral([a, b, c, d], new_dih, incr=False, axis_pos="bond", threshold=0.1)`	Rotates the component on `c`'s side of bond `b–c` around that central bond; default pivot is its midpoint
`axis_pos="term"`	Rotates atoms beyond `d` about a line through `d` parallel to `b–c`; this legacy mode leaves the selected quartet's dihedral unchanged
`threshold`	Nonnegative angular tolerance in degrees; skips changes at or below this tolerance without shifting planar targets away from 0° or ±180°

The distance and angle setters require exactly two or three Python integer positions respectively. set_dihedral accepts four distinct Python or NumPy integer positions within the state. The setters return None and do not update graph topology, cached fragments, or Atom.x/y/z. set_dihedral requires finite target and tolerance values, validates axis_pos, and uses a finite symmetric conn_mat indexed by atom position. It preserves connectivity while selecting the atoms to move and rejects attempts to rotate a central bond in a ring before editing coordinates or connectivity.

The distance and angle setters still use the legacy determine_carried_atoms helper. That helper temporarily removes the chosen bond from conn_mat, selects the component containing its second atom, and normally restores the bond. If the graph remains connected, it prints Matrix still connected, returns None, and currently fails to restore the removed bond. Use those setters and the helper only on bonds whose removal separates the relevant components, and keep a copy. set_dihedral does not use this helper; it does not support ring-closure editing either.

Maintain atoms and connectivity explicitly¶

State.add_atoms(atom) appends one atom and State.remove_atoms(atom) removes one object. Neither method updates counts, IDs, or graphs. After structural changes, explicitly synchronize these fields before calling the older geometry tools:

from matmalib import Atom, State, StateEditor
import networkx as nx

state = State(id="pair", NAtoms=1, atoms=[
    Atom(coord=[0.0, 0.0, 0.0], type="C", index=0)
])
added = Atom(coord=[1.4, 0.0, 0.0], type="C", index=1)
state.add_atoms(added)
state.NAtoms = len(state.atoms)
for index, atom in enumerate(state.atoms):
    atom.id = index
state.connectivity_matrix()

StateEditor.disconnect_atoms(state, 0, 1)
print("Disconnected:", state.conn_mat[0, 1])
StateEditor.connect_atoms(state, 0, 1)
print("Connected:", state.conn_mat[0, 1])
# These helpers change only conn_mat. Rebuild graph when needed.
state.graph = nx.from_numpy_array(state.conn_mat)
state.Nmols = nx.number_connected_components(state.graph)
state.system_split()
print("Fragments:", len(state.mols))

state.remove_atoms(added)
state.NAtoms = len(state.atoms)
for index, atom in enumerate(state.atoms):
    atom.id = index
state.connectivity_matrix()
state.system_split()
print("Remaining atoms:", state.NAtoms)

For custom bonds, initialize a symmetric (NAtoms, NAtoms) matrix, keep its diagonal zero, and create graph with nx.from_numpy_array. connect_atoms and disconnect_atoms mutate the two symmetric entries only; adjacent_atoms returns the nonzero column indices in a row.

Inspect a structure in a notebook¶

State.show(width=600, height=600) creates a Py3Dmol stick rendering, zooms to the structure, and displays it. It requires an environment that renders Py3Dmol's HTML/JavaScript. It returns None rather than a viewer handle.

from matmalib import Space

state = Space.from_xyz("tests/data/glucose.xyz").states[0]
state.show(width=700, height=500)

The browser tutorial runner executes Python and Matplotlib examples; this notebook viewer example should be run in Jupyter. Continue to reading and writing data or working with ensembles.