Preparing aromatic participants#

Separating covalent cycles from aromatic participant recognition.

Prepare complete ring memberships before a geometric analysis. This recipe does not calculate pi–pi interactions. Both tools are experimental.

Added in version 1.0.0.

1. Declare a synthetic ring graph#

Supply atom and bond flags explicitly. No hydrogens, geometry, protonation, or force-field parameters are inferred; this is a declared graph fixture.

import molsysmt as msm
import pandas as pd
from molsysmt.native import Topology

# A synthetic ring graph with explicitly supplied chemical flags.
molsys = Topology(n_atoms=6)
molsys.atoms['atom_type'] = ['C'] * 6
molsys.bonds = pd.DataFrame({
    'atom1_index': [0, 1, 2, 3, 4, 5],
    'atom2_index': [1, 2, 3, 4, 5, 0],
    'bond_type': ['covalent'] * 6,
    'is_aromatic': [True] * 6,
})
msm.set(molsys, element='atom', atom_is_aromatic=[True] * 6)

2. Prepare chemical participants#

Connectivity perception and chemical recognition return sparse atom-index memberships. They report the selected chemical state and the actual software versions used.

cycles = msm.topology.get_rings(molsys, assume_complete_connectivity=True)
aromatic = msm.physchem.get_aromatic_rings(molsys, assume_complete_connectivity=True)
assert cycles['atom_indices'].tolist() == aromatic['atom_indices'].tolist()
assert aromatic['atom_offsets'].tolist() == [0, 6]
states = msm.convert(molsys, to_form='molsysmt.ChemicalStatesDict')
assert msm.physchem.get_aromatic_rings(
    states, assume_complete_connectivity=True)['atom_indices'].tolist() == aromatic['atom_indices'].tolist()
aromatic['definition'], aromatic['rule_version'], aromatic['software']

('stored_aromatic_bonds',
 'stored_aromatic_bond_cycles@1',
 {'molsysmt': '0.21.0+606.ga03eb4bf6', 'networkx': '3.6.1'})

3. Select complete memberships#

Use a list of original atom indices including all atoms of each intersected ring. Fused rings can share atoms, so choosing only one subring can cut its neighbor. Choose complete fused components when necessary.

selected = msm.physchem.get_aromatic_rings(
    molsys, selection=list(range(6)), assume_complete_connectivity=True)
participants = [selected['atom_indices'][a:b] for a, b in zip(
    selected['atom_offsets'][:-1], selected['atom_offsets'][1:])]
participants

[array([0, 1, 2, 3, 4, 5])]

4. Fit geometric planes#

Attach explicitly unitful coordinates to a MolSys carrying the declared chemistry, then pass complete memberships to the general structure tool. This synthetic hexagon is a geometric control, not a molecular simulation. Neither planarity nor its normal supplies aromaticity; the chemical flags supplied above do that.

import numpy as np
from molsysmt import pyunitwizard as puw

molsys = msm.convert(molsys, to_form='molsysmt.MolSys')
angles = np.arange(6) * np.pi / 3
xyz = np.zeros((1, 6, 3))
xyz[0, :, :2] = .14 * np.column_stack((np.cos(angles), np.sin(angles)))
msm.append_structures(molsys, from_molecular_system=puw.quantity(xyz, 'nm'))
planes = msm.structure.get_least_squares_plane(molsys, selection=participants)
print('Plane normal:', planes['normals'][0, 0])
print('RMS planarity:', planes['rms_deviation'])

Plane normal: [0. 0. 1.]
RMS planarity: [[0.0]] nanometer

5. Preserve declarations through H5MSM#

Save chemical metadata in the existing H5MSM 0.5 domains. The reader recomputes this small declared-bond basis; the participant dictionary is not itself an attached interaction analysis.

from tempfile import TemporaryDirectory
from pathlib import Path

with TemporaryDirectory() as directory:
    filename = str(Path(directory) / 'rings.h5msm')
    msm.convert(molsys, to_form=filename)
    restored = msm.physchem.get_aromatic_rings(
        filename, assume_complete_connectivity=True)
assert restored['atom_indices'].tolist() == aromatic['atom_indices'].tolist()
assert restored['chemical_state_index'] == aromatic['chemical_state_index']

See also