Saving hydrogen-bond analyses#

Calculating, naming, persisting and querying attributed hydrogen bonds.

Compare definitions without overwriting results or losing the calculation provenance. This controlled water pair provides a small, reproducible pipeline with evaluated-empty structures, named analyses and explicit source indices.

Added in version 1.0.0.

Calculating named definitions#

Each method produces an independent analysis. Attach it explicitly after calculation; the MolSys domain validates compatible axes without authenticating independent origins.

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

molsys = msm.convert(Chem.AddHs(Chem.MolFromSmiles('O.O')), to_form='molsysmt.MolSys')
xyz = np.array([[0, 0, 0], [.28, 0, 0], [.1, 0, 0], [0, .1, 0], [.28, .1, 0], [.28, 0, .1]])
coordinates = np.repeat(xyz[None], 3, axis=0)
coordinates[1, [1, 4, 5]] += [0, .4, 0]
molsys.structures.append(coordinates=puw.quantity(coordinates, 'nm'))

analyses = {method: msm.interactions.hbonds.get_hbonds(molsys, method=method, pbc=False)
            for method in ('baker_hubbard', 'cpptraj', 'prolif')}
molsys.interactions = {**molsys.interactions, **analyses}
print({name: result.n_interactions for name, result in molsys.interactions.items()})

{'baker_hubbard': 2, 'cpptraj': 2, 'prolif': 2}

Persisting the system and provenance#

Public conversion writes H5MSM 0.5. It retains method references, original producer versions, measurement units, coverage and stable occurrence indices for this analysis version. The installed reader version does not replace the saved producer version.

import tempfile
from pathlib import Path
workspace = tempfile.TemporaryDirectory()
path = str(Path(workspace.name) / 'hydrogen_bonds.h5msm')
msm.convert(molsys, to_form=path)
restored = msm.convert(path, to_form='molsysmt.MolSys')
for name, result in analyses.items():
    loaded = restored.interactions[name]
    assert loaded.software == result.software
    assert loaded.parameters == result.parameters
    assert loaded.to_dict()['occurrence_indices'].tolist() == result.to_dict()['occurrence_indices'].tolist()
print(list(restored.interactions))

['baker_hubbard', 'cpptraj', 'prolif']

Querying and remapping selected structures#

Query the visible frame without recalculating. A frame evaluated without interactions differs from one never evaluated. Extraction preserves requested order and remaps local indices while retaining source maps. Removing a participant removes its relation.

analysis = restored.interactions['baker_hubbard']
print(analysis.query(structure_indices=[2, 0, 2], atom_indices=[2]).to_dict()['occurrence_indices'])
assert analysis.query(structure_indices=1).n_interactions == 0
subset = msm.extract(restored, structure_indices=[2, 0])
print(subset.interactions['baker_hubbard'].structure_source_indices.tolist())
without_acceptor = msm.extract(restored, selection=[0, 2, 3])
assert without_acceptor.interactions['baker_hubbard'].n_interactions == 0

[1 0]
[2, 0]

Calculating directly from file blocks#

Index selections read chemical information once and only projected coordinate blocks. Saved analyses are not required for a fresh calculation; accepted output remains in RAM.

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

2 [0, 2]