Molecule Graph¶

The Atomistic class is the fundamental graph representation in MolPy. Unlike the Frame class which is optimized for bulk array operations, simulations, Atomistic is optimized for chemistry.

It treats molecules as mathematical graphs where: Nodes are Atoms, Atom. and Edges are Bonds, Bond..

This tutorial covers Graph Theory Basics: Understanding the data structure., CRUD Operations: Creating, Reading, Updating, and Deleting atoms/bonds., Traversal: Efficiently querying connectivity., Topology: Auto-generating angles and dihedrals., and Subgraphs: Extracting chemical components.

1. Initialization¶

An Atomistic object starts as an empty container. It holds: atoms: A collection, Entities of Atom objects. and bonds: A collection, Entities of Bond objects..

In [1]:

import molpy as mp

# Initialize an empty molecule
# We can pass arbitrary metadata, stored in the object's properties
mol = mp.Atomistic(metadata={"name": "ethanol"})

print(f"Molecule created. Atoms: {len(mol.atoms)}, Bonds: {len(mol.bonds)}")

import molpy as mp # Initialize an empty molecule # We can pass arbitrary metadata, stored in the object's properties mol = mp.Atomistic(metadata={"name": "ethanol"}) print(f"Molecule created. Atoms: {len(mol.atoms)}, Bonds: {len(mol.bonds)}")

Molecule created. Atoms: 0, Bonds: 0

2. Atoms¶

Creating Atoms¶

Use def_atom() to add a node. You can pass any arbitrary key-value pairs, **kwargs. Common attributes include: name: Unique string ID, optional but recommended for debugging., element: Chemical symbol., x, y, z: Coordinates, Angstroms., charge: Partial charge, elementary charge units., and type: Force field atom type..

In [2]:

# Define the carbon backbone
c1 = mol.def_atom(element="C", name="C1", x=0.0, y=0.0, z=0.0, charge=-0.2)
c2 = mol.def_atom(element="C", name="C2", x=1.54, y=0.0, z=0.0, charge=-0.1)

# Define the oxygen
o1 = mol.def_atom(element="O", name="O1", x=2.0, y=1.4, z=0.0, charge=-0.6)

# Define the hydroxyl hydrogen
h1 = mol.def_atom(element="H", name="HO", x=2.9, y=1.4, z=0.0, charge=0.4)

# Accessing atoms by index or iteration
print(f"Added atoms: {[a.get('name') for a in mol.atoms]}")

# Define the carbon backbone c1 = mol.def_atom(element="C", name="C1", x=0.0, y=0.0, z=0.0, charge=-0.2) c2 = mol.def_atom(element="C", name="C2", x=1.54, y=0.0, z=0.0, charge=-0.1) # Define the oxygen o1 = mol.def_atom(element="O", name="O1", x=2.0, y=1.4, z=0.0, charge=-0.6) # Define the hydroxyl hydrogen h1 = mol.def_atom(element="H", name="HO", x=2.9, y=1.4, z=0.0, charge=0.4) # Accessing atoms by index or iteration print(f"Added atoms: {[a.get('name') for a in mol.atoms]}")

Added atoms: ['C1', 'C2', 'O1', 'HO']

Reading and Updating Atoms¶

Atoms behave like Python dictionaries. You can access or modify properties using bracket notation [] or .get().

Note: The Atom object is a reference. Modifying it updates the graph immediately.

In [3]:

# Accessing data
print(f"C1 Charge: {c1['charge']}")

# Updating data
c1['type'] = 'opls_135'  # Assigning an OPLS atom type
print(f"C1 Type: {c1.get('type')}")

# Bulk update example
for atom in mol.atoms:
    if atom.get('element') == 'C':
        atom['hybridization'] = 'sp3'

print(f"C2 Hybridization: {c2.get('hybridization')}")

# Accessing data print(f"C1 Charge: {c1['charge']}") # Updating data c1['type'] = 'opls_135' # Assigning an OPLS atom type print(f"C1 Type: {c1.get('type')}") # Bulk update example for atom in mol.atoms: if atom.get('element') == 'C': atom['hybridization'] = 'sp3' print(f"C2 Hybridization: {c2.get('hybridization')}")

C1 Charge: -0.2
C1 Type: opls_135
C2 Hybridization: sp3

3. Bonds¶

Creating Bonds¶

Use def_bond, atom1, atom2 to connect nodes. Like atoms, bonds can store attributes, e.g., order, type, k, r0.

In [4]:

# Connect the heavy atoms
b1 = mol.def_bond(c1, c2, type="C-C", order=1)
b2 = mol.def_bond(c2, o1, type="C-O", order=1)
b3 = mol.def_bond(o1, h1, type="O-H", order=1)

print(f"Created {len(mol.bonds)} bonds.")
# The Bond object exposes endpoints via 'itom' and 'jtom' properties
print(f"Bond 1: {b1.itom.get('name')} -- {b1.jtom.get('name')} (Type: {b1.get('type')})")

# Connect the heavy atoms b1 = mol.def_bond(c1, c2, type="C-C", order=1) b2 = mol.def_bond(c2, o1, type="C-O", order=1) b3 = mol.def_bond(o1, h1, type="O-H", order=1) print(f"Created {len(mol.bonds)} bonds.") # The Bond object exposes endpoints via 'itom' and 'jtom' properties print(f"Bond 1: {b1.itom.get('name')} -- {b1.jtom.get('name')} (Type: {b1.get('type')})")

Created 3 bonds.
Bond 1: C1 -- C2 (Type: C-C)

Checking Connectivity¶

To check connectivity or get neighbors, use the Atomistic container methods, e.g. get_neighbors. The atoms themselves do not store connectivity, the container does.

In [5]:

# Check if C2 is connected to O1
neighbors_c2 = mol.get_neighbors(c2)

is_connected = o1 in neighbors_c2
print(f"Is C2 connected to O1? {is_connected}")

# Check if C2 is connected to O1 neighbors_c2 = mol.get_neighbors(c2) is_connected = o1 in neighbors_c2 print(f"Is C2 connected to O1? {is_connected}")

Is C2 connected to O1? True

4. Graph Traversal¶

We can traverse the graph using mol.get_neighbors, atom.

In [6]:

# Who is attached to C2?
print(f"Neighbors of {c2.get('name')}:")
for neighbor in mol.get_neighbors(c2):
    print(f"  -> {neighbor.get('name')} ({neighbor.get('element')})")

# What bonds is C2 involved in?
# We can iterate over all bonds in the molecule and check endpoints
print(f"\nBonds of {c2.get('name')}:")
for bond in mol.bonds:
    if c2 in bond.endpoints:
        # Find the partner atom
        partner = bond.itom if bond.jtom is c2 else bond.jtom
        print(f"  -> {bond.get('type')} bond to {partner.get('name')}")

# Who is attached to C2? print(f"Neighbors of {c2.get('name')}:") for neighbor in mol.get_neighbors(c2): print(f" -> {neighbor.get('name')} ({neighbor.get('element')})") # What bonds is C2 involved in? # We can iterate over all bonds in the molecule and check endpoints print(f"\nBonds of {c2.get('name')}:") for bond in mol.bonds: if c2 in bond.endpoints: # Find the partner atom partner = bond.itom if bond.jtom is c2 else bond.jtom print(f" -> {bond.get('type')} bond to {partner.get('name')}")

Neighbors of C2:
  -> C1 (C)
  -> O1 (O)

Bonds of C2:
  -> C-C bond to C1
  -> C-O bond to O1

Advanced Query: Finding Neighbor Depth¶

We can easily write recursive functions to find atoms within N hops using mol.get_neighbors().

In [7]:

def get_n_hop_neighbors(struct, start_atom, n=2):
    visited = {start_atom}
    current_shell = {start_atom}
    
    for _ in range(n):
        next_shell = set()
        for atom in current_shell:
            # Use struct.get_neighbors
            for neighbor in struct.get_neighbors(atom):
                if neighbor not in visited:
                    visited.add(neighbor)
                    next_shell.add(neighbor)
        current_shell = next_shell
    
    return visited - {start_atom}

neighborhood = get_n_hop_neighbors(mol, c1, n=2)
# Use set comprehension for display
names = {a.get('name') for a in neighborhood}
print(f"Atoms within 2 hops of C1: {names}")

def get_n_hop_neighbors(struct, start_atom, n=2): visited = {start_atom} current_shell = {start_atom} for _ in range(n): next_shell = set() for atom in current_shell: # Use struct.get_neighbors for neighbor in struct.get_neighbors(atom): if neighbor not in visited: visited.add(neighbor) next_shell.add(neighbor) current_shell = next_shell return visited - {start_atom} neighborhood = get_n_hop_neighbors(mol, c1, n=2) # Use set comprehension for display names = {a.get('name') for a in neighborhood} print(f"Atoms within 2 hops of C1: {names}")

Atoms within 2 hops of C1: {'O1', 'C2'}

5. Topology Generation¶

Molecular Dynamics requires not just bonds, but also Angles, 3-body interactions and Dihedrals, 4-body interactions.

MolPy can automatically derive these from the connectivity graph using get_topo().

In [8]:

# Initially, no angles/dihedrals exist
print(f"Before: {len(mol.angles)} angles, {len(mol.dihedrals)} dihedrals")

# Generate topology
mol.get_topo(gen_angle=True, gen_dihe=True)

print(f"After:  {len(mol.angles)} angles, {len(mol.dihedrals)} dihedrals")

# Initially, no angles/dihedrals exist print(f"Before: {len(mol.angles)} angles, {len(mol.dihedrals)} dihedrals") # Generate topology mol.get_topo(gen_angle=True, gen_dihe=True) print(f"After: {len(mol.angles)} angles, {len(mol.dihedrals)} dihedrals")

Before: 0 angles, 0 dihedrals
After:  2 angles, 1 dihedrals

Inspecting Topology¶

Each Angle/Dihedral object references the atoms involved.

In [9]:

print("Angles found:")
for ang in mol.angles:
    # angles have itom, jtom, ktom properties or we can use endpoints
    names = '-'.join([a.get('name') for a in ang.endpoints])
    print(f"  {names}")

print("\nDihedrals found:")
for dihe in mol.dihedrals:
    names = '-'.join([a.get('name') for a in dihe.endpoints])
    print(f"  {names}")

print("Angles found:") for ang in mol.angles: # angles have itom, jtom, ktom properties or we can use endpoints names = '-'.join([a.get('name') for a in ang.endpoints]) print(f" {names}") print("\nDihedrals found:") for dihe in mol.dihedrals: names = '-'.join([a.get('name') for a in dihe.endpoints]) print(f" {names}")

Angles found:
  C1-C2-O1
  C2-O1-HO

Dihedrals found:
  C1-C2-O1-HO

6. Subgraphs and Modification¶

Removing Atoms¶

Removing an atom automatically removes all connected bonds, angles, and dihedrals to maintain graph consistency.

In [10]:

# Let's remove the H1 atom
print(f"Removing {h1.get('name')}...")
mol.remove_entity(h1)

print(f"Atom count: {len(mol.atoms)}")
print(f"Bond count: {len(mol.bonds)}")

# Verify bonds are gone
neighbors_o1 = [n.get('name') for n in mol.get_neighbors(o1)]
print(f"Neighbors of O1: {neighbors_o1}")

# Let's remove the H1 atom print(f"Removing {h1.get('name')}...") mol.remove_entity(h1) print(f"Atom count: {len(mol.atoms)}") print(f"Bond count: {len(mol.bonds)}") # Verify bonds are gone neighbors_o1 = [n.get('name') for n in mol.get_neighbors(o1)] print(f"Neighbors of O1: {neighbors_o1}")

Removing HO...
Atom count: 3
Bond count: 2
Neighbors of O1: ['C2', 'C1', 'C2']

Copying Molecules¶

You often need to clone a molecule, e.g., to create multiple solvent molecules. Use .copy() to get a deep copy with new unique IDs.

In [11]:

mol_copy = mol.copy()

# Verify independence
# Note: We need to access atoms in the copy. They are new objects.
# Let's find C1 in the copy (it will have the same name)
c1_copy = [a for a in mol_copy.atoms if a.get('name') == 'C1'][0]
c1_copy['name'] = "C1_modified"

print(f"Original C1 name: {c1.get('name')}")
print(f"Copy C1 name:     {c1_copy.get('name')}")

mol_copy = mol.copy() # Verify independence # Note: We need to access atoms in the copy. They are new objects. # Let's find C1 in the copy (it will have the same name) c1_copy = [a for a in mol_copy.atoms if a.get('name') == 'C1'][0] c1_copy['name'] = "C1_modified" print(f"Original C1 name: {c1.get('name')}") print(f"Copy C1 name: {c1_copy.get('name')}")

Original C1 name: C1
Copy C1 name:     C1_modified

Summary¶

Atomistic is the core editable molecular graph in MolPy. It provides a unified representation of atoms and bonds and supports direct graph-level modifications, including defining atoms and bonds and removing existing entities. Atoms and bonds behave as dictionary-like containers, allowing arbitrary properties to be attached and queried without imposing a fixed schema.

Connectivity within the molecular graph is accessed through neighbor traversal, enabling local and global graph operations directly from atomic relationships. Based on the bond graph, Atomistic can automatically derive higher-order topology such as angles and dihedrals via get_topo(), ensuring that the structure remains consistent and ready for downstream molecular dynamics workflows without manual bookkeeping.

To see how this graph-based model scales to large systems, continue with User Guide: Polymer Stepwise, which demonstrates how thousands of atoms are assembled automatically into polymer structures.