molify

https://github.com/zincware/molify
Admin
molify is a molecular structure interface that integrates RDKit, ASE, Packmol, and NetworkX to...
Tokens:21,190
Snippets:118
Trust Score:7.4
License:Apache-2.0
Update:10 months ago
Show doc for...
Tokens:
Raw
### Install molify from PyPI

Source: https://github.com/zincware/molify/blob/main/docs/source/index.rst

Installs the molify package using pip from the Python Package Index. This is the standard method for installing Python libraries. Ensure you have pip and a virtual environment set up.

```console
(.venv) $ pip install molify
```

--------------------------------

### Install molify from Source using UV

Source: https://github.com/zincware/molify/blob/main/docs/source/index.rst

Installs molify from source using the uv package manager, including recursive submodules. This method is recommended for development. Requires git and uv to be installed.

```console
(.venv) $ git clone --recursive https://github.com/zincware/molify
(.venv) $ cd molify
(.venv) $ uv sync
(.venv) $ source .venv/bin/activate
```

--------------------------------

### Convert ASE Molecules to RDKit for Visualization in Python

Source: https://github.com/zincware/molify/blob/main/docs/source/ase_tools.ipynb

This example shows how to convert molecules from the ASE G2 database into RDKit molecules for visualization and analysis. It iterates through the first five molecules in the database and displays them. Requires 'molify', 'ase', and 'rdkit'.

```python
from ase.collections import g2
import molify
from IPython.display import display

# Show a few molecules from the G2 database
for idx, (atoms, name) in enumerate(zip(g2, g2.names)):
    if idx >= 5:  # Show first 5
        break

    print(f"\n{name}:")
    # Note: These molecules don't have connectivity, so we need suggestions
    # We'll learn more about this in the NetworkX tools section
    mol = molify.ase2rdkit(atoms, suggestions=[])
    display(mol)
```

--------------------------------

### Multiple matches in propanol using SMILES and SMARTS

Source: https://github.com/zincware/molify/blob/main/docs/source/atom_selection.ipynb

Shows how `match_substructure` handles multiple occurrences of a pattern within a molecule, using propanol as an example. It finds all carbon atoms and a specific C-C-O fragment, returning all matches.

```python
# Create propanol with multiple carbon atoms
propanol = molify.smiles2atoms("CCCO")
propanol_mol = molify.ase2rdkit(propanol, suggestions=["CCCO"])

# Find all carbon atoms
carbon_matches = molify.match_substructure(propanol_mol, "[#6]")
print(f"All carbon matches in propanol: {carbon_matches}")
print(f"Total carbon atoms found: {len(carbon_matches)}")

# Find the C-C-O pattern
cco_pattern = molify.match_substructure(propanol_mol, "CCO")
print(f"\nC-C-O pattern matches: {cco_pattern}")

# Visualize
Draw.MolToImage(propanol_mol, size=(300, 200))
```

--------------------------------

### Calculate Cubic Box Dimensions for System Setup

Source: https://context7.com/zincware/molify/llms.txt

Computes the dimensions of a cubic simulation box required to achieve a specific target density for a given set of molecules. This utility function, calculate_box_dimensions, is crucial for setting up simulations with appropriate system sizes. It requires molecular objects and the target density as input.

```python
from molify import smiles2conformers
from molify.utils import calculate_box_dimensions
import ase.units

# Create molecules
water = smiles2conformers("O", 1)
ethanol = smiles2conformers("CCO", 1)

# Example usage (assuming calculate_box_dimensions is available and used)
# box_dimensions = calculate_box_dimensions([water, ethanol], density=1000)
# print(f"Calculated box dimensions: {box_dimensions}")

```

--------------------------------

### Import molify and RDKit modules

Source: https://github.com/zincware/molify/blob/main/docs/source/atom_selection.ipynb

Initializes the molify library and imports necessary RDKit modules for chemical structure manipulation and drawing.

```python
from rdkit.Chem import CombineMols, Draw

import molify
```

--------------------------------

### Initialize Git Submodules (Console)

Source: https://github.com/zincware/molify/blob/main/docs/source/index.rst

This command updates and initializes Git submodules. It's particularly useful if the `--recursive` flag was not used during the initial `git clone` operation. This ensures all necessary submodules are fetched and prepared for use within the project. No specific inputs are required beyond the command itself.

```console
(.venv) $ git submodule update --init --recursive
```

--------------------------------

### Create a Mixture of Water and Ethanol with molify

Source: https://github.com/zincware/molify/blob/main/docs/source/packmol_tools.ipynb

Illustrates how to create a mixed molecular system containing both water and ethanol. It demonstrates generating conformers for multiple molecule types and packing them into a single box with a specified density and composition.

```python
# Generate conformers for both molecule types
water = molify.smiles2conformers("O", numConfs=10)
ethanol = molify.smiles2conformers("CCO", numConfs=10)

print(f"Water conformers: {len(water)}")
print(f"Ethanol conformers: {len(ethanol)}")

# Create mixture: 5 water + 5 ethanol
mixture = molify.pack(
    data=[water, ethanol],
    counts=[5, 5],
    density=800,  # kg/m³
)

print(f"Mixture: {mixture}")
print(f"Cell: {mixture.cell.lengths()} Å")
print(f"Total atoms: {len(mixture)}")
```

--------------------------------

### Import molify and RDKit

Source: https://github.com/zincware/molify/blob/main/docs/source/packmol_tools.ipynb

Initializes the necessary libraries for using molify and RDKit functionalities. RDKit is often used for molecular manipulation and visualization, and molify leverages it for conformer generation.

```python
from rdkit.Chem import Draw

import molify
```

--------------------------------

### Import Molify and NetworkX Libraries

Source: https://github.com/zincware/molify/blob/main/docs/source/networkx_tools.ipynb

Imports necessary libraries for molecular manipulation and graph analysis. Includes ASE for atomic structures, NetworkX for graph operations, IPython for display, RDKit for chemical drawing, and Molify for core functionalities.

```python
import ase.build
import networkx as nx
from IPython.display import display
from rdkit.Chem import Draw

import molify
```

--------------------------------

### Create a Simple Water Box with molify

Source: https://github.com/zincware/molify/blob/main/docs/source/packmol_tools.ipynb

Demonstrates the basic usage of molify's `pack()` function to create a periodic box filled with water molecules at a specified density. It shows conformer generation and box creation, including outputting box properties.

```python
# Generate water conformers
water = molify.smiles2conformers("O", numConfs=10)

print(f"Generated {len(water)} water conformers")
print(f"Each conformer: {water[0]}")

# Pack 10 water molecules into a box
water_box = molify.pack(
    data=[water],  # List of conformer lists
    counts=[10],  # Number of each molecule type
    density=1000,  # Density in kg/m³
)

print(f"Water box: {water_box}")
print(f"Cell dimensions: {water_box.cell.lengths()} Å")
print(f"Total atoms: {len(water_box)}")
print(f"Volume: {water_box.get_volume():.2f} Å³")
```

--------------------------------

### Create and Analyze Benzene Molecule using RDKit and NetworkX

Source: https://github.com/zincware/molify/blob/main/docs/source/rdkit_tools.ipynb

Demonstrates creating a benzene molecule from SMILES notation using RDKit, converting it to a NetworkX graph using molify, and then identifying cycles and bond orders within the aromatic ring. Requires RDKit, NetworkX, and molify libraries.

```python
# Create benzene
benzene_mol = Chem.MolFromSmiles("c1ccccc1")
benzene_mol = Chem.AddHs(benzene_mol)

print("Benzene molecule:")
display(benzene_mol)

# Convert to graph
benzene_graph = molify.rdkit2networkx(benzene_mol)
print(
    f"\nNodes: {benzene_graph.number_of_nodes()}, "
    f"Edges: {benzene_graph.number_of_edges()}"
)
```

```python
# Find all cycles in benzene
cycles = nx.cycle_basis(benzene_graph)
print(f"Number of cycles found: {len(cycles)}\n")

# The aromatic ring should be the largest cycle
aromatic_ring = max(cycles, key=len)
print(f"Aromatic ring (6-membered): {aromatic_ring}")

# Check bond orders in the aromatic ring
print("\nBond orders in the aromatic ring:")
for i in range(len(aromatic_ring)):
    node1 = aromatic_ring[i]
    node2 = aromatic_ring[(i + 1) % len(aromatic_ring)]
    if benzene_graph.has_edge(node1, node2):
        bond_order = benzene_graph[node1][node2]["bond_order"]
        print(f"  C{node1} - C{node2}: {bond_order} (aromatic)")
```

--------------------------------

### Specify and Handle Bond Orders with NetworkX and RDKit

Source: https://context7.com/zincware/molify/llms.txt

Demonstrates how to explicitly set bond orders when creating molecular graphs using NetworkX and converting them to RDKit molecules. It also shows how to handle unknown bond orders by providing suggestions for determination.

```python
from molify import networkx2rdkit
import networkx as nx
from rdkit import Chem

# Specify bond orders explicitly
graph = nx.Graph()
graph.add_edge(0, 1, bond_order=1.0)
graph.add_edge(1, 2, bond_order=2.0)

mol = networkx2rdkit(graph)
print(f"SMILES: {Chem.MolToSmiles(mol)}")
print(f"Bonds: {[(b.GetBeginAtomIdx(), b.GetEndAtomIdx(), b.GetBondTypeAsDouble()) for b in mol.GetBonds()]}')

# Handle unknown bond orders with suggestions
graph_unknown = nx.Graph()
graph_unknown.add_node(0, atomic_number=6, charge=0,
                       position=[0.0, 0.0, 0.0])
graph_unknown.add_node(1, atomic_number=8, charge=0,
                       position=[1.2, 0.0, 0.0])
graph_unknown.add_edge(0, 1, bond_order=None)

# Provide hints for bond order determination
mol_determined = networkx2rdkit(graph_unknown, suggestions=["C=O"])
print(f"Determined structure: {Chem.MolToSmiles(mol_determined)}")
```

--------------------------------

### Import molify and ASE Libraries

Source: https://github.com/zincware/molify/blob/main/docs/source/ase_tools.ipynb

Imports necessary modules from ASE, IPython.display, and molify for working with molecular structures and simulations. This is a foundational step for using molify's functionalities.

```python
import ase.build
from IPython.display import display

import molify
```

--------------------------------

### Visualize Molecular Mixture using RDKit

Source: https://github.com/zincware/molify/blob/main/docs/source/packmol_tools.ipynb

Converts the packed molecular system (ASE Atoms object) into an RDKit molecule object for visualization. This allows for graphical representation and inspection of the generated mixture.

```python
# Convert to RDKit for visualization
mixture_mol = molify.ase2rdkit(mixture)

# Visualize
Draw.MolToImage(mixture_mol, size=(600, 400))
```

--------------------------------

### Calculate System Density using Molify Utilities (Python)

Source: https://context7.com/zincware/molify/llms.txt

Computes the density of an ASE Atoms object using the 'calculate_density' utility from the molify library. This function works with any ASE Atoms object that has periodic boundary conditions (PBC) enabled. It demonstrates creating a packed system and then calculating its density, along with verifying the relationship between mass, volume, and density.

```python
from molify import smiles2conformers, pack
from molify.utils import calculate_density

# Create packed system
water = smiles2conformers("O", 1)
box = pack([water], [100], density=1000, seed=42)

# Calculate density
computed_density = calculate_density(box)

print(f"Target density: 1000 kg/m³")
print(f"Computed density: {computed_density:.0f} kg/m³")

# Works with any ASE Atoms with PBC
import ase
import numpy as np
test_atoms = ase.Atoms(
    'H2O',
    positions=[[0, 0, 0], [1, 0, 0], [0, 1, 0]],
    cell=[10, 10, 10],
    pbc=True
)
test_density = calculate_density(test_atoms)
print(f"Test system density: {test_density:.2f} kg/m³")

# Relationship: density = mass / volume
print(f"Total mass: {sum(box.get_masses()):.2f} amu")
print(f"Volume: {box.get_volume():.2f} Å³")
```

--------------------------------

### Build Molecular Systems with Packmol

Source: https://github.com/zincware/molify/blob/main/docs/source/index.rst

Uses Packmol to construct molecular systems, such as solvated molecules or mixtures, at a specified density. This function is valuable for preparing input for molecular dynamics simulations. Requires molify and Packmol.

```python
import molify

# Generate conformers for water and ethanol
water = molify.smiles2conformers("O", numConfs=10)
etoh = molify.smiles2conformers("CCO", numConfs=10)

# Pack 5 water and 5 ethanol molecules into a box at 800 kg/m³
box = molify.pack(
   data=[water, etoh], counts=[5, 5], density=800
)
```

--------------------------------

### Create and Visualize Toluene Molecule with Highlights

Source: https://context7.com/zincware/molify/llms.txt

Generates a toluene molecule from SMILES notation, identifies aromatic and aliphatic carbons, and visualizes the molecule with different highlights for these carbon types. It uses molify functions like smiles2atoms, ase2rdkit, match_substructure, and visualize_selected_molecules. The output is a PIL Image object.

```python
from molify import smiles2atoms, ase2rdkit, match_substructure, visualize_selected_molecules

# Create toluene molecule
toluene = smiles2atoms("Cc1ccccc1")
mol = ase2rdkit(toluene, suggestions=["Cc1ccccc1"])

# Find aromatic and aliphatic carbons
aromatic = match_substructure(mol, "c", hydrogens="exclude")
aliphatic = match_substructure(mol, "[C;!c]", hydrogens="exclude")

# Flatten tuples to lists
aromatic_flat = [idx for match in aromatic for idx in match]
aliphatic_flat = [idx for match in aliphatic for idx in match]

# Visualize with different highlights (different colors)
img = visualize_selected_molecules(
    mol,
    aromatic_flat,
    aliphatic_flat,
    mols_per_row=2,
    sub_img_size=(300, 300),
    legends=["Toluene with highlights"],
    alpha=0.6
)

# Display returns PIL Image object
# img.show()  # Uncomment to display
print(f"Generated image: {type(img)}")

# Visualize without highlights
simple_img = visualize_selected_molecules(mol)
```

--------------------------------

### Analyze Molecular Cycles using NetworkX

Source: https://github.com/zincware/molify/blob/main/docs/source/networkx_tools.ipynb

Illustrates finding cycles within a molecule's graph representation using NetworkX. It converts an ASE structure to a NetworkX graph, prints basic graph statistics, and then uses `nx.cycle_basis` to identify and visualize all cycles.

```python
# Create serotonin (has cycles)
serotonin = molify.smiles2atoms("C1=CC2=C(C=C1O)C(=CN2)CCN")
serotonin_mol = molify.ase2rdkit(serotonin)

display(serotonin_mol)

# Analyze with NetworkX
graph = molify.ase2networkx(serotonin)
print(f"\nSerotonin: {graph.number_of_nodes()} atoms, {graph.number_of_edges()} bonds")
molify.draw_molecular_graph(graph)

# Find all cycles
cycles = nx.cycle_basis(graph)
print(f"Found {len(cycles)} cycles:\n")

for idx, cycle in enumerate(cycles):
    print(f"Cycle {idx + 1} ({len(cycle)} atoms): {cycle}")
    # Visualize each cycle
    display(Draw.MolToImage(serotonin_mol, highlightAtoms=cycle, size=(300, 200)))
```

--------------------------------

### Pack Molecules into Periodic Box using molify and Packmol

Source: https://context7.com/zincware/molify/llms.txt

Builds multi-component molecular systems within a periodic box using Packmol integration. Users can specify the number of molecules for each component, desired density, and simulation parameters like seed, tolerance, and box shape (cubic or non-cubic). The output is an ASE Atoms object representing the packed system.

```python
from molify import pack, smiles2conformers

# Prepare molecular components with multiple conformers
water = smiles2conformers("O", numConfs=2, randomSeed=1)
ethanol = smiles2conformers("CCO", numConfs=5, randomSeed=2)

# Pack 50 water + 20 ethanol molecules at 1000 kg/m³
density = 1000  # kg/m³
packed_system = pack(
    data=[water, ethanol],
    counts=[50, 20],
    density=density,
    seed=42,
    tolerance=2.0,
    pbc=True
)

print(packed_system)
# Output: Atoms(symbols='C40H180O70', pbc=True, cell=[...])

print(f"Box dimensions: {packed_system.cell.diagonal()}")
print(f"Total atoms: {len(packed_system)}")
print(f"Volume: {packed_system.get_volume():.2f} Å³")
print(f"Density: {sum(packed_system.get_masses()) * 1.66054e-27 / (packed_system.get_volume() * 1e-30):.0f} kg/m³")

# Create non-cubic box with 2:1:1 aspect ratio
rectangular_box = pack(
    data=[water],
    counts=[100],
    density=1000,
    ratio=(2.0, 1.0, 1.0)
)
print(f"Non-cubic cell: {rectangular_box.cell.diagonal()}")
```

--------------------------------

### Analyze Molecules in a Mixture with NetworkX

Source: https://github.com/zincware/molify/blob/main/docs/source/networkx_tools.ipynb

Shows how to analyze a mixture of molecules (water and ethanol) using NetworkX. It creates a packed system, converts it to a graph, and then uses `nx.connected_components` to identify individual molecules and analyze their sizes.

```python
# Create a water/ethanol mixture
water = molify.smiles2conformers("O", numConfs=10) ethanol = molify.smiles2conformers("CCO", numConfs=10)
box = molify.pack(data=[water, ethanol], counts=[5, 5], density=800)

# Convert to graph
box_graph = molify.ase2networkx(box)
print(
    f"System: {box_graph.number_of_nodes()} atoms, {box_graph.number_of_edges()} bonds"
)

from collections import Counter

# Find connected components (individual molecules)
components = list(nx.connected_components(box_graph))
print(f"Found {len(components)} separate molecules")

# Analyze molecule sizes
sizes = Counter(len(comp) for comp in components)
print("\nMolecule sizes:")
for size, count in sorted(sizes.items()):
    molecule_type = "water (H2O)" if size == 3 else "ethanol (C2H6O)"
    print(f"  {count} molecules with {size} atoms ({molecule_type})")
```

--------------------------------

### Build Periodic Boxes with Packmol

Source: https://github.com/zincware/molify/blob/main/README.md

Uses the Packmol interface to construct periodic boxes of molecules at a specified density. This function takes a list of molecular conformers and their desired counts, then uses Packmol to arrange them in a simulation box.

```python
from molify import pack, smiles2conformers

water = smiles2conformers("O", 2) ethanol = smiles2conformers("CCO", 5)

density = 1000  # kg/m^3
box = pack([water, ethanol], [7, 5], density)
print(box)
>>> Atoms(symbols='C10H44O12', pbc=True, cell=[8.4, 8.4, 8.4])
```

--------------------------------

### Calculate Density of Packed System

Source: https://github.com/zincware/molify/blob/main/docs/source/packmol_tools.ipynb

Demonstrates how to calculate the actual density of a generated molecular box using the `calculate_density` utility function from `molify.utils`. This function takes an ASE Atoms object as input.

```python
# Calculate actual density
from molify.utils import calculate_density

actual_density = calculate_density(mixture)
print(f"Density: {actual_density:.1f} kg/m³")
```

--------------------------------

### Create Ammonia Molecule with ASE

Source: https://github.com/zincware/molify/blob/main/docs/source/networkx_tools.ipynb

Demonstrates creating an ammonia molecule using the `ase.build.molecule` function. It also checks if the resulting ASE Atoms object contains connectivity information, which is often missing for simple inorganic molecules built this way.

```python
# Create ammonia from ASE (no connectivity info)
ammonia = ase.build.molecule("NH3")
print(f"Ammonia: {ammonia}")
print(f"Has connectivity: {'connectivity' in ammonia.info}")
```

--------------------------------

### Mapped SMILES for Precise Atom Selection and Ordering

Source: https://github.com/zincware/molify/blob/main/docs/source/atom_selection.ipynb

Illustrates the use of mapped SMILES patterns with `mapped_only=True` in `molify.match_substructure` to select specific atoms and ensure they are returned in ascending order of their map numbers. This provides deterministic atom ordering for various downstream applications.

```python
# Using mapped SMILES to select only specific atoms
# This pattern matches the carbon-carbon-oxygen chain,
# but with mapped_only=True, only returns the mapped carbons
mapped_pattern = "[C:1][C:2]O"
mapped_indices = molify.match_substructure(
    ethanol_mol, mapped_pattern, mapped_only=True
)
print(f"Mapped carbon indices: {mapped_indices}")

# Compare with mapped_only=False which returns all matched atoms
unmapped_indices = molify.match_substructure(
    ethanol_mol, mapped_pattern, mapped_only=False
)
print(f"All matched atoms: {unmapped_indices}")

# Demonstrate atom ordering control with alanine dipeptide
aladip = molify.smiles2atoms("CC(=O)NC(C)C(=O)NC")
aladip_mol = molify.ase2rdkit(aladip, suggestions=["CC(=O)NC(C)C(=O)NC"])

# Select atoms in map order 1, 2, 3, 4
indices_1234 = molify.match_substructure(
    aladip_mol, "CC(=O)N[C:1]([C:2])[C:3](=O)[N:4]C", mapped_only=True
)
print(f"\nMap order 1,2,3,4: {indices_1234}")

# Select the same atoms but in different map order 4, 3, 2, 1
indices_4321 = molify.match_substructure(
    aladip_mol, "CC(=O)N[C:4]([C:3])[C:2](=O)[N:1]C", mapped_only=True
)
print(f"Map order 4,3,2,1: {indices_4321}")
```

--------------------------------

### Visualize Molecular Highlights and Custom Coloring

Source: https://context7.com/zincware/molify/llms.txt

Introduces the `visualize_selected_molecules` function for generating 2D molecular visualizations. This function allows highlighting specific atoms or bonds and applying custom coloring schemes.

```python
from molify import smiles2atoms, ase2rdkit, match_substructure
from molify import visualize_selected_molecules

# Placeholder for visualization example - actual code depends on visualization backend.
```

--------------------------------

### Round-trip Conversion: ASE to RDKit

Source: https://github.com/zincware/molify/blob/main/docs/source/rdkit_tools.ipynb

Demonstrates a lossless conversion process by converting an ASE Atoms object (representing aspirin) back into an RDKit Mol object using molify. This verifies that chemical information is preserved during the conversion. Requires molify and RDKit.

```python
# ASE → RDKit (round-trip)
aspirin_roundtrip = molify.ase2rdkit(aspirin_atoms)

print("Original molecule:")
display(aspirin_mol)

print("\nRound-trip molecule (ASE → RDKit):")
display(aspirin_roundtrip)
```

--------------------------------

### Visualizing Atom Selections

Source: https://github.com/zincware/molify/blob/main/docs/source/atom_selection.ipynb

Introduces the `visualize_selected_molecules` function for highlighting selected atoms within molecular structures. It supports visualizing multiple selection sets with distinct colors and requires flattening match results from `match_substructure` into a simple list of indices.

```python
# Create toluene molecule for visualization examples
toluene_smiles = "Cc1ccccc1"
toluene_mol = molify.ase2rdkit(
    molify.smiles2atoms(toluene_smiles), suggestions=["Cc1ccccc1"]
)

# Select different atom types
aromatic_carbons = molify.match_substructure(toluene_mol, "c", hydrogens="exclude")
methyl_carbon = molify.match_substructure(toluene_mol, "[C;!c]", hydrogens="exclude")

print(f"Aromatic carbons: {aromatic_carbons}")
print(f"Methyl carbon: {methyl_carbon}")

# Flatten matches for visualization
aromatic_flat = [idx for match in aromatic_carbons for idx in match]
methyl_flat = [idx for match in methyl_carbon for idx in match]

```

--------------------------------

### Convert ASE Atoms to NetworkX Graph

Source: https://context7.com/zincware/molify/llms.txt

Generates NetworkX Graph representations from ASE Atoms objects. Connectivity is automatically detected based on covalent radii and a scaling factor, or uses explicit connectivity if provided in `atoms.info`. Handles periodic systems.

```python
from molify import smiles2atoms, ase2networkx

# Create molecular structure
atoms = smiles2atoms("CCO")

# Convert to graph with automatic connectivity detection
graph = ase2networkx(atoms, pbc=False, scale=1.2)

print(f"Nodes: {len(graph.nodes)}")
print(f"Edges: {len(graph.edges)}")

# Access spatial information
for node, attrs in graph.nodes(data=True):
    print(f"Atom {node}: element={attrs['atomic_number']}, "
          f"position={attrs['position']}")

# Check if explicit connectivity was used
if 'connectivity' in atoms.info:
    print("Used explicit connectivity from atoms.info")
    bond_orders = [e[2]['bond_order'] for e in graph.edges(data=True)]
    print(f"Bond orders: {bond_orders}")
else:
    print("Used distance-based connectivity detection")

# For periodic systems
from molify import pack, smiles2conformers
water_conf = smiles2conformers("O", 1)
periodic_box = pack([water_conf], [10], density=1000)
periodic_graph = ase2networkx(periodic_box, pbc=True)
print(f"Periodic system: {len(periodic_graph.nodes)} atoms")
```

--------------------------------

### Match carbon and hydroxyl oxygen atoms in ethanol

Source: https://github.com/zincware/molify/blob/main/docs/source/atom_selection.ipynb

Demonstrates basic substructure matching using SMARTS patterns to find all carbon atoms and the hydroxyl oxygen atom in an ethanol molecule. It returns the atom indices for each match.

```python
# Match carbon atoms using SMARTS pattern
carbon_matches = molify.match_substructure(ethanol_mol, "[#6]")
print(f"Carbon atom matches: {carbon_matches}")
print(f"Number of carbon matches: {len(carbon_matches)}")

# Match the hydroxyl oxygen
hydroxyl_matches = molify.match_substructure(ethanol_mol, "[#8]")
print(f"\nHydroxyl oxygen matches: {hydroxyl_matches}")
```

--------------------------------

### Demonstrate Scale Parameter for Bond Detection in Python

Source: https://github.com/zincware/molify/blob/main/docs/source/ase_tools.ipynb

This snippet demonstrates how the `scale` parameter affects bond detection in Molify's `ase2networkx()` function. It iterates through different scale values and prints the number of bonds detected for each. Requires the 'molify' library and ASE.

```python
import molify
from ase.build import molecule

ammonia = molecule('NH3')

# Try different scale values
scales = [0.7, 1.2, 3]

for scale in scales:
    graph = molify.ase2networkx(ammonia, scale=scale)
    _ = molify.draw_molecular_graph(graph)
    print(f"scale={scale}: {graph.number_of_edges()} bonds detected")
```

--------------------------------

### Group Substructure Matches by Molecular Fragment

Source: https://context7.com/zincware/molify/llms.txt

Demonstrates how to use `group_matches_by_fragment` to organize substructure matches based on disconnected molecular components in a multi-component system. This is useful for analyzing complex mixtures.

```python
from rdkit import Chem
from molify import match_substructure, group_matches_by_fragment

# Create molecule with multiple fragments
mol = Chem.MolFromSmiles("CCO.CF.CCCC")  # Ethanol + fluoromethane + butane
mol = Chem.AddHs(mol)

# Find all carbons
matches = match_substructure(mol, "[#6]")
print(f"Total carbon matches: {matches}")

# Group by fragment
grouped = group_matches_by_fragment(mol, matches)
print(f"Carbons per fragment: {[len(g) for g in grouped]}")
# Output: [2, 1, 4] for three fragments

for i, group in enumerate(grouped):
    print(f"Fragment {i}: carbon indices {group}")

# Each fragment's carbons are separate
assert len(grouped) == 3  # Three disconnected molecules
assert sum(len(g) for g in grouped) == len(matches)  # All matches accounted for
```

--------------------------------

### Compress Molecular System to Higher Density

Source: https://context7.com/zincware/molify/llms.txt

Demonstrates compressing a molecular system to a target density while optionally preserving molecular geometry. It utilizes molify functions smiles2conformers, pack, and compress. The input is a molecular system, and the output is a compressed system with updated cell dimensions and volume.

```python
from molify import smiles2conformers, pack, compress

# Create initial system
water = smiles2conformers("O", 1)
initial_box = pack([water], [50], density=800, seed=42)

print(f"Initial density: {800} kg/m³")
print(f"Initial cell: {initial_box.cell.diagonal()}")
print(f"Initial volume: {initial_box.get_volume():.2f} Å³")

# Compress to higher density
compressed = compress(initial_box, density=1000, freeze_molecules=False)

print(f"Target density: {1000} kg/m³")
print(f"Compressed cell: {compressed.cell.diagonal()}")
print(f"Compressed volume: {compressed.get_volume():.2f} Å³")

# Verify volume ratio matches density ratio
ratio = initial_box.get_volume() / compressed.get_volume()
print(f"Volume ratio: {ratio:.2f} (expected: {1000/800:.2f})")

# Compress with frozen molecular geometry
frozen_compressed = compress(
    initial_box,
    density=1000,
    freeze_molecules=True
)
print("Frozen compression: molecular geometries preserved")
# Requires connectivity information in atoms.info
```

--------------------------------

### Convert RDKit Mol and ASE Atoms to NetworkX Graph

Source: https://github.com/zincware/molify/blob/main/docs/source/index.rst

Generates a NetworkX Graph representation from either an RDKit Mol object or an ASE Atoms object. This enables network analysis on molecular structures. Requires RDKit, ASE, and molify.

```python
import molify
import ase.build
from rdkit import Chem

etoh_mol = Chem.MolFromSmiles("CCO")
nhi_atoms = ase.build.molecule("NH3")

etoh_graph = molify.rdkit2networkx(etoh_mol)
nhi_graph = molify.ase2networkx(nh3_atoms)
```

--------------------------------

### Handling Conformer Reuse in molify

Source: https://github.com/zincware/molify/blob/main/docs/source/packmol_tools.ipynb

Explains and demonstrates the behavior of molify's `pack()` function when the requested number of molecules (`counts`) exceeds the number of generated conformers. In such cases, conformers are reused.

```python
water = molify.smiles2conformers("O", numConfs=5)
box = molify.pack(
    data=[water],
    counts=[10],  # Requesting 10 water molecules will use each conformer twice
    density=997 # kg/m^3
)
```

--------------------------------

### Iterate Over Molecular Fragments in Multi-Component Systems

Source: https://context7.com/zincware/molify/llms.txt

Shows how to decompose a system containing multiple molecules into its individual components using `iter_fragments`. This is applicable to systems created by packing or simple concatenation.

```python
from molify import smiles2atoms, pack, smiles2conformers, iter_fragments
import ase

# Create multi-component system
water = smiles2conformers("O", 1)
ethanol = smiles2conformers("CCO", 1)
mixed_system = pack([water, ethanol], [3, 2], density=1000, seed=42)

# Iterate over fragments
fragments = list(iter_fragments(mixed_system))
print(f"Found {len(fragments)} molecules")

for i, frag in enumerate(fragments):
    print(f"Molecule {i}: {frag.get_chemical_formula()}")

# Count molecule types
formulas = [f.get_chemical_formula() for f in fragments]
print(f"Water: {formulas.count('H2O')}")
print(f"Ethanol: {formulas.count('C2H6O')}")

# Works with simple concatenation too
combined = smiles2atoms("CCO") + smiles2atoms("O")
combined_frags = list(iter_fragments(combined))
assert len(combined_frags) == 2
```

--------------------------------

### Convert ASE Atoms to NetworkX Graph and Back

Source: https://github.com/zincware/molify/blob/main/docs/source/networkx_tools.ipynb

Demonstrates converting an ASE Atoms object (representing a molecule) into a NetworkX graph and then back into an ASE Atoms object. This process preserves positional and chemical information and verifies connectivity preservation.

```python
# Create a simple molecule and convert to graph
water = molify.smiles2atoms("O")
water_graph = molify.ase2networkx(water)

print(
    f"Graph: {water_graph.number_of_nodes()} nodes, "
    f"{water_graph.number_of_edges()} edges"
)

# Convert back to ASE
water_back = molify.networkx2ase(water_graph)
print(f"\nBack to ASE: {water_back}")
print(f"Connectivity preserved: {'connectivity' in water_back.info}")
print(f"Number of bonds: {len(water_back.info['connectivity'])}")
```

--------------------------------

### Specific SMARTS matching for sp3 carbons and hydroxyls

Source: https://github.com/zincware/molify/blob/main/docs/source/atom_selection.ipynb

Illustrates using more specific SMARTS patterns to match sp3 hybridized carbon atoms and oxygen atoms with exactly one bonded hydrogen (hydroxyl group) in an ethanol molecule.

```python
# Use SMARTS for more specific pattern matching
# Match only sp3 carbon atoms (aliphatic carbons with single bonds)
sp3_carbons = molify.match_substructure(ethanol_mol, "[C;X4]")
print(f"sp3 carbon matches: {sp3_carbons}")

# Match oxygen with exactly one hydrogen (hydroxyl group)
hydroxyl_oxygen = molify.match_substructure(ethanol_mol, "[OH1]")
print(f"Hydroxyl oxygen: {hydroxyl_oxygen}")
```

--------------------------------

### Visualize Highlighted Molecules with Molify

Source: https://github.com/zincware/molify/blob/main/docs/source/atom_selection.ipynb

This snippet visualizes a molecule (e.g., toluene) and highlights specific selected atoms, such as aromatic or methyl carbons. It takes a molecule object and lists of atom indices to highlight as input, returning an image object.

```python
import molify

# Assume toluene_mol, aromatic_flat, and methyl_flat are defined
# For example:
toluene_smiles = "c1ccccc1C"
toluene_mol = molify.ase2rdkit(molify.smiles2atoms(toluene_smiles))
aromatic_flat = [1, 2, 3, 4, 5, 6] # Example indices for aromatic carbons
methyl_flat = [0] # Example index for methyl carbon

img = molify.visualize_selected_molecules(toluene_mol, aromatic_flat, methyl_flat)
print(img)
```

--------------------------------

### Create and visualize ethanol molecule

Source: https://github.com/zincware/molify/blob/main/docs/source/atom_selection.ipynb

Creates an ethanol molecule from SMILES string, prints its atom count and chemical formula, and converts it to an RDKit Mol object for visualization and further processing. This conversion step is crucial for substructure matching.

```python
# Create an ethanol molecule for basic examples
ethanol = molify.smiles2atoms("CCO")
print(f"Ethanol has {len(ethanol)} atoms")
print(f"Chemical formula: {ethanol.get_chemical_formula()}")

# Convert to RDKit Mol for substructure matching
ethanol_mol = molify.ase2rdkit(ethanol, suggestions=["CCO"])

# Visualize the structure
Draw.MolToImage(ethanol_mol, size=(300, 200))
```

--------------------------------

### Match Molecular Substructures with Flexible Hydrogen Handling

Source: https://context7.com/zincware/molify/llms.txt

Uses the `match_substructure` function to find occurrences of molecular patterns defined by SMARTS. Supports flexible hydrogen handling ('include', 'isolated') and atom mapping for precise selection.

```python
from molify import smiles2atoms, ase2rdkit, match_substructure

# Create molecule
atoms = smiles2atoms("CCCO")  # Propanol
mol = ase2rdkit(atoms, suggestions=["CCCO"])

# Find all carbon atoms
carbons = match_substructure(mol, "[#6]")
print(f"Carbon atoms: {carbons}")
# Output: ((0,), (1,), (2,))

# Find C-O bonds including hydrogens
co_with_h = match_substructure(mol, "CO", hydrogens="include")
print(f"C-O with hydrogens: {co_with_h}")
# Returns indices of C, O, and all H atoms bonded to them

# Get only hydrogens on matched atoms
h_only = match_substructure(mol, "CO", hydrogens="isolated")
print(f"Hydrogens on C-O: {h_only}")

# Use atom mapping to select specific atoms
mapped = match_substructure(mol, "[C:1][C:2]O", mapped_only=True)
print(f"Mapped carbons in order: {mapped}")
# Returns in map order: C:1, C:2

# Find aromatic carbons in benzene
benzene = ase2rdkit(smiles2atoms("c1ccccc1"))
aromatic = match_substructure(benzene, "c")
print(f"Found {len(aromatic)} aromatic carbons")
```

--------------------------------

### Convert ASE to NetworkX

Source: https://github.com/zincware/molify/blob/main/docs/source/networkx_tools.ipynb

Converts an ASE Atoms object to a NetworkX graph. This conversion may result in bond orders being set to None.

```python
ammonia_graph = molify.ase2networkx(ammonia)

print("\nGraph edge attributes:")
for u, v, data in ammonia_graph.edges(data=True):
    print(f"  N({u})-H({v}): bond_order = {data['bond_order']}")
```

--------------------------------

### Convert NetworkX Graph to RDKit Molecule with Bond Orders

Source: https://github.com/zincware/molify/blob/main/docs/source/networkx_tools.ipynb

Explains and demonstrates the conversion of a NetworkX graph, which may contain explicit bond order information, into an RDKit molecule object. This is crucial when the graph originates from sources like RDKit itself.

```python
# Create molecule with known bonds
acetone = molify.smiles2atoms("CC(=O)C")  # Has a C=O double bond
acetone_graph = molify.ase2networkx(acetone)

# Check bond orders
print("Bond orders in graph:")
for u, v, data in list(acetone_graph.edges(data=True))[:5]:
    print(f"  {u}-{v}: {data['bond_order']}")

# Convert to RDKit
acetone_mol = molify.networkx2rdkit(acetone_graph)
print("\nConversion successful:")
display(acetone_mol)
```

--------------------------------

### Convert ASE Atoms to RDKit Molecule

Source: https://context7.com/zincware/molify/llms.txt

Transforms ASE Atoms objects back into RDKit Molecules. It can infer bond orders from existing connectivity or molecular geometry. Supports automatic aromaticity and double bond detection, and accepts suggestions for complex cases.

```python
from molify import smiles2atoms, ase2rdkit
from rdkit import Chem

# Create ASE structure
atoms = smiles2atoms("CCO")

# Convert back to RDKit
mol = ase2rdkit(atoms)

print(f"Number of atoms: {mol.GetNumAtoms()}")
print(f"Number of bonds: {mol.GetNumBonds()}")
print(f"SMILES: {Chem.MolToSmiles(mol)}")
# Output: CCO

# Check bond orders are preserved
for bond in mol.GetBonds():
    print(f"Bond {bond.GetBeginAtomIdx()}-{bond.GetEndAtomIdx()}: "
          f"order={bond.GetBondTypeAsDouble()}")

# Use suggestions for complex molecules with ambiguous bonding
benzene = smiles2atoms("c1ccccc1")
mol_benzene = ase2rdkit(benzene, suggestions=["c1ccccc1"])
print(f"Aromatic bonds: {sum(1 for b in mol_benzene.GetBonds() if b.GetIsAromatic())}")
```

--------------------------------

### Convert NetworkX to RDKit with Bond Order Determination

Source: https://github.com/zincware/molify/blob/main/docs/source/networkx_tools.ipynb

Converts a NetworkX graph to an RDKit molecule. This process automatically triggers bond order determination if not already present. The 'suggestions' parameter can be used to provide known SMILES strings to aid the determination.

```python
# Convert to RDKit - triggers automatic bond order determination
ammonia_mol = molify.networkx2rdkit(ammonia_graph, suggestions=[])

display(ammonia_mol)
```

```python
SMILES = "OP(=O)(O)OP(=O)(O)OC=CC1C(O)C(OC2C(C(O)C(N3C=NC4=C3N=CN=C4N)O2)O)O1"

molecule = molify.smiles2atoms(SMILES)
molecule.info.pop("connectivity", None)  # Remove connectivity info
graph = molify.ase2networkx(molecule)

# the automatic bond order determination can fail for complex molecules
try:
    mol = molify.networkx2rdkit(graph, suggestions=[])
except ValueError as e:
    print(f"Conversion failed: {e}")

# but if we know the molecule, we can provide it as a suggestion
# You can also provide a list of suggestions for multiple molecules or multiple guesses
mol = molify.networkx2rdkit(graph, suggestions=[SMILES])
display(mol)
# For very large molecules, this can take a while!
```

--------------------------------

### Convert ASE Atoms to RDKit Mol

Source: https://github.com/zincware/molify/blob/main/docs/source/index.rst

Converts an ASE Atoms object to an RDKit Mol object. This allows RDKit's cheminformatics tools to be used with ASE-represented structures. Dependencies include ASE and molify.

```python
import molify
import ase.build

nhi = ase.build.molecule("NH3")
mol = molify.ase2rdkit(nh3)
```

--------------------------------

### Extract Substructures as ASE Atoms Objects

Source: https://context7.com/zincware/molify/llms.txt

Explains how to use `get_substructures` to extract specific molecular patterns or fragments from a molecule and return them as ASE Atoms objects. This allows for further analysis and manipulation of substructures.

```python
from molify import smiles2atoms, get_substructures

# Create propanol molecule
propanol = smiles2atoms("CCCO")

# Extract all carbon-containing fragments
carbons = get_substructures(propanol, "[#6]", suggestions=["CCCO"])

print(f"Extracted {len(carbons)} carbon atoms")
for i, carbon in enumerate(carbons):
    print(f"Carbon {i}: {carbon.get_chemical_formula()}, "
          f"position: {carbon.positions[0]}")

# Extract specific functional groups
alcohols = get_substructures(propanol, "CO", suggestions=["CCCO"])
print(f"Alcohol groups: {len(alcohols)}")

# Extract larger substructures
ethyl_groups = get_substructures(propanol, "CC", suggestions=["CCCO"])
print(f"Ethyl fragments: {len(ethyl_groups)}")
```

--------------------------------

### Analyze and Visualize Aspirin Functional Groups with Molify

Source: https://github.com/zincware/molify/blob/main/docs/source/atom_selection.ipynb

This Python snippet demonstrates how to identify and visualize multiple functional groups (acetyl ester, carbonyl oxygens, carboxylic acid hydrogen) within aspirin using Molify's substructure matching capabilities. It converts SMILES to a molecule object, performs matches, flattens the results, and then visualizes the molecule with distinct highlights and legends.

```python
# Aspirin (acetylsalicylic acid) contains multiple functional groups
aspirin_smiles = "CC(=O)Oc1ccccc1C(=O)O"
aspirin_mol = molify.ase2rdkit(
    molify.smiles2atoms(aspirin_smiles), suggestions=[aspirin_smiles]
)

# Select different functional groups
# Acetyl ester group with hydrogens (using mapped_only to get specific atoms)
acetyl = molify.match_substructure(
    aspirin_mol,
    "[C:1][C:2](=O)[O:4]c1ccccc1C(=O)O",
    hydrogens="include",
    mapped_only=True,
)

# Carbonyl oxygens
cabornyls = molify.match_substructure(aspirin_mol, "[O:1]=C", mapped_only=True)

# Carboxylic acid hydrogens
acid_h = molify.match_substructure(aspirin_mol, "CO", hydrogens="isolated")

# Flatten for visualization
acetyl_flat = [idx for match in acetyl for idx in match]
cabornyls_flat = [idx for match in cabornyls for idx in match]
acid_h_flat = [idx for match in acid_h for idx in match]

# Visualize with multiple functional group highlights
img_result = molify.visualize_selected_molecules(
    aspirin_mol,
    acetyl_flat,
    cabornyls_flat,
    acid_h_flat,
    mols_per_row=1,
    legends=["Aspirin functional group analysis"],
)
print(img_result)
```

--------------------------------

### Verify Connectivity Preservation in Packed Box

Source: https://github.com/zincware/molify/blob/main/docs/source/packmol_tools.ipynb

Verifies that the connectivity information, including bond orders and atom connectivity, is accurately preserved after packing molecules into a box using molify. This is crucial for maintaining molecular integrity in simulations.

```python
# Connectivity is preserved!
print(f"Has connectivity: {'connectivity' in water_box.info}")
print(f"Number of bonds: {len(water_box.info['connectivity'])}")
print(f"Expected bonds: {10 * 2} (10 molecules x 2 O-H bonds each)")
```
Tokens:
Raw