Saving aromatic interaction analyses#

Calculating, naming, persisting and querying sparse aromatic ring observations.

Keep a geometric analysis beside the system that supplied its indices. This controlled two-ring ensemble has parallel, absent and edge-to-face observations. The cutoffs define an experimental geometric criterion, not an energy model. The SMILES conversion requires the optional RDKit dependency.

Added in version 1.0.0.

Preparing the chemical system#

Create declared aromatic chemistry and a controlled coordinate ensemble.

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

# RDKit supplies explicit aromatic atom/bond flags through this form conversion.
topology = msm.convert('smiles:c1ccccc1.c1ccccc1', to_form='molsysmt.Topology')
molsys = msm.convert(topology, to_form='molsysmt.MolSys')
phase = np.arange(6) * np.pi / 3
ring = .14 * np.column_stack((np.cos(phase), np.sin(phase), np.zeros(6)))
coordinates = np.array([
    np.concatenate((ring, ring + [0, 0, .35])),
    np.concatenate((ring, ring + [0, 0, 2.])),
    np.concatenate((ring, ring[:, [2, 1, 0]] + [0, 0, .4])),
])
molsys.structures.append(coordinates=puw.quantity(coordinates, 'nm'))

def calculate(source=molsys, **options):
    settings = dict(distance_threshold='.6 nm', angle_threshold='30 degrees',
                    offset_threshold='.2 nm', planarity_threshold='.02 nm', pbc=False)
    settings.update(options)
    return msm.interactions.pi_pi.get_pi_pi_interactions(source, **settings)

analysis = calculate()
print(analysis.n_interactions, analysis.occurrence_structures.tolist())
print(analysis.measure_units)
assert analysis.n_interactions == 2
assert not molsys.interactions

2 [0, 2]
{'distance': 'nm', 'plane_angle': 'radians', 'offset_a': 'nm', 'offset_b': 'nm', 'rms_deviation_a': 'nm', 'rms_deviation_b': 'nm', 'max_deviation_a': 'nm', 'max_deviation_b': 'nm'}

Combining separate domains#

A pipeline may hold chemistry and coordinates separately. Their pair declares matching atom indices. The same detector consumes that composite, and public conversion builds a topology-free destination for the named result. Several chemical states would still need a declared per-structure assignment when using chemical_state='structure'.

chemical_states = msm.convert(molsys.chemical_states,
                              to_form='molsysmt.ChemicalStatesDict')
domains = [chemical_states, molsys.structures]
composite_analysis = calculate(domains)
np.testing.assert_array_equal(composite_analysis.occurrence_structures,
                              analysis.occurrence_structures)
np.testing.assert_allclose(composite_analysis.measurements['distance'],
                           analysis.measurements['distance'], atol=1e-12, rtol=0)
molsys = msm.convert(domains, to_form='molsysmt.MolSys')
analysis = composite_analysis
assert molsys.topology is None
print('Topology-free system:', molsys.structures.n_atoms,
      'atoms;', molsys.structures.n_structures, 'structures')

Topology-free system: 12 atoms; 3 structures

Saving the named analysis#

Attachment declares correspondence with the system atom and structure axes. H5MSM 0.5 preserves both coordinates and the named sparse analysis.

from pathlib import Path
from tempfile import TemporaryDirectory

molsys.interactions = {**molsys.interactions, 'aromatic_contacts': analysis}
with TemporaryDirectory() as directory:
    filename = str(Path(directory) / 'aromatic.h5msm')
    msm.convert(molsys, to_form=filename)
    loaded = msm.convert(filename, to_form='molsysmt.MolSys')
    saved = loaded.interactions['aromatic_contacts']
    with msm.configure.context(chunk_size=1):
        from_file = calculate(source=filename, heavy_mode='force', structure_indices=[2, 0, 2, 1])
    np.testing.assert_array_equal(saved.query().to_dict()['occurrence_indices'],
                                  analysis.query().to_dict()['occurrence_indices'])
    np.testing.assert_array_equal(from_file.occurrence_structures, analysis.occurrence_structures)
    for name in analysis.measurements:
        np.testing.assert_allclose(saved.measurements[name], analysis.measurements[name])
        np.testing.assert_allclose(from_file.measurements[name], analysis.measurements[name])
    assert saved.parameters == analysis.parameters
    assert saved.software == analysis.software
print(saved.software)
print(saved.evaluated_structure_indices.tolist())

{'molsysmt': '0.22.4+24.g42b487869', 'networkx': '3.7'}
[0, 1, 2]

Querying the visible structure#

An evaluated empty frame is not an unevaluated frame. Queries preserve occurrence indices for selection and inspection by a client such as MolSysViewer.

visible = saved.query(structure_indices=[2], atom_indices=[0], mode='involving_selection').to_dict()
empty = saved.query(structure_indices=[1]).to_dict()
assert len(visible['occurrence_indices']) == 1
assert empty['evaluated_structure_indices'].tolist() == [1]
assert len(empty['occurrence_indices']) == 0
print('Visible occurrence indices:', visible['occurrence_indices'])

Visible occurrence indices: [1]

Saving a standalone layer#

A standalone analysis needs an explicit declaration of correspondence before attachment to another system; matching axis sizes alone does not prove identity.

with TemporaryDirectory() as directory:
    filename = str(Path(directory) / 'analysis_only.h5msm')
    msm.h5msm.write_layers(filename, interactions={'aromatic_contacts': analysis})
    layers = msm.h5msm.read_layers(filename, layers='interactions')
    standalone = layers['interactions']['aromatic_contacts']
assert standalone.software == analysis.software
assert standalone.n_interactions == analysis.n_interactions
print('Sparse numeric bytes:', standalone.numeric_nbytes)

Sparse numeric bytes: 680