Saving cation-pi analyses#

Calculating, naming, persisting and querying attributed and proposed cation–π methods.

Keep the method and its scientific provenance attached to the molecular system. This controlled benzene/ammonium ensemble includes an evaluated empty frame. The workflow uses only temporary files and does not install the reference ProLIF package.

Added in version 1.0.0.

Calculating separately named methods#

The default reproduces ProLIF 2.2.2’s CationPi definition. The custom rule uses four explicit cutoffs and remains a proposal under comparison, not a demonstrated scientific improvement. Attach both explicitly without overwriting other analyses.

import numpy as np
import molsysmt as msm
from molsysmt import pyunitwizard as puw
from rdkit import Chem

molsys = msm.convert(Chem.MolFromSmiles('c1ccccc1.[NH4+]'), to_form='molsysmt.MolSys')
phase = np.arange(6) * np.pi / 3
ring = np.column_stack((.14*np.cos(phase), .14*np.sin(phase), np.zeros(6)))
points = [[0, 0, .35], [.25, 0, .35], [0, 0, 2], [0, 0, -.35], [.1, 0, .35]]
coordinates = np.asarray([np.vstack((ring, point)) for point in points])
molsys.structures.append(coordinates=puw.quantity(coordinates, 'nm'))
known = msm.interactions.cation_pi.get_cation_pi_interactions(molsys, pbc=False)
custom = msm.interactions.cation_pi.get_cation_pi_interactions(molsys, method='centroid_angle_offset', distance_threshold='.45 nm', angle_threshold='30 degrees', offset_threshold='.2 nm', planarity_threshold='.02 nm', pbc=False)
molsys.interactions = {**molsys.interactions, 'prolif': known, 'proposal': custom}
print({name: result.n_interactions for name, result in molsys.interactions.items()})

{'prolif': 3, 'proposal': 3}

Saving chemistry, structures and named observations#

Public conversion writes H5MSM 0.5. Method-reference version and producer software are different facts; both survive loading. Declared bracket hydrogen counts and SMARTS ring order also survive. Query one visible frame without recalculation.

import tempfile
from pathlib import Path
workspace = tempfile.TemporaryDirectory()
path = str(Path(workspace.name) / 'system.h5msm')
msm.convert(molsys, to_form=path)
restored = msm.convert(path, to_form='molsysmt.MolSys')
print(list(restored.interactions))
print(restored.interactions['prolif'].parameters['method_reference']['version'])
print(restored.interactions['prolif'].software)
print(restored.interactions['prolif'].query(structure_indices=2).n_interactions)
assert restored.interactions['prolif'].parameters == known.parameters
assert msm.get(restored, element='atom', n_explicit_hydrogens=True) == msm.get(molsys, element='atom', n_explicit_hydrogens=True)

['prolif', 'proposal']
2.2.2
{'molsysmt': '0.21.0+606.ga03eb4bf6', 'rdkit': '2025.09.5'}
0

Computing directly from projected file blocks#

Atom-index selections or ‘all’ permit the bounded H5MSM route. It reads chemistry once and only required coordinate blocks, without saved analyses. The accepted analysis remains in RAM; smaller coordinate blocks do not guarantee smaller process RSS or a faster run. Rich file selections use bounded eager source materialization.

with msm.configure.context(chunk_size=1):
    projected = msm.interactions.cation_pi.get_cation_pi_interactions(path, structure_indices=[4, 0, 2, 4], heavy_mode='force', pbc=False)
print(projected.evaluated_structure_indices.tolist(), projected.occurrence_structures.tolist())
print(projected.execution_records[0]['details']['execution_chunks'])

[0, 2, 4] [0, 4]
3

Extracting frames and remapping observations#

Source maps retain original atom/structure indices. Named analyses accompany an extraction; they are not recomputed. Removing one atom of a compound ring discards that relation explicitly. remap can reorder axes within the declared source but does not embed a subsystem into a larger destination.

subset = msm.extract(restored, structure_indices=[4, 0, 2])
print(subset.interactions['prolif'].structure_source_indices.tolist())
print(subset.interactions['prolif'].query(atom_indices=[6]).n_interactions)
reordered = known.remap(atom_indices=[6, 5, 4, 3, 2, 1, 0], structure_indices=[4, 0])
print(reordered.atom_source_indices.tolist(), reordered.query(atom_indices=[0]).n_interactions)

[4, 0, 2]
2
[6, 5, 4, 3, 2, 1, 0] 2

Writing an independent interactions layer#

A typed dictionary supports in-memory interchange. H5MSM stores typed columns, not a Python object per occurrence. A file containing only analyses has no coordinates. Attaching it to another system requires declaring axis correspondence; matching axis lengths alone does not establish atom identity.

columns = msm.convert(known, to_form='molsysmt.InteractionsDict')
assert msm.convert(columns, to_form='molsysmt.Interactions').n_interactions == known.n_interactions
analysis_path = str(Path(workspace.name) / 'analysis.h5msm')
msm.h5msm.write_layers(analysis_path, interactions={'prolif': known})
loaded = msm.h5msm.read_layers(analysis_path, layers='interactions')['interactions']['prolif']
assert loaded.parameters == known.parameters
assert loaded.software == known.software
print(loaded.query(structure_indices=[0, 4], atom_indices=[6]).to_dict()['occurrence_indices'])
workspace.cleanup()

[0 2]