Adding hydrogens to a fixed chemical state#

Materializing a prepared hydrogen inventory while preserving the source pose.

Use molsysmt.build.add_missing_hydrogens() with mode='fixed_chemical_state', pH=None and engine='RDKit' when you have already chosen the component’s chemistry. This experimental mode returns a new native MolSys, regardless of the supported input form. Set return_report=True to receive {'molecular_system': ..., 'report': ...}. RDKit is optional and loaded only when this engine is requested. An unavailable engine raises; no fallback runs.

API documentation

Added in version 1.0.0.

Choosing chemistry before geometry#

Supply one isolated connected component, one chemical state and one finite coordinate structure. Explicitly extract other states/structures first; selecting a structure alone does not authorize their removal. Required assignments are element symbols (atom_type), complete covalent connectivity, supported bond orders, atom/bond aromaticity, formal charges, zero unpaired electrons and stored implicit/explicit H counts. This slice supports H, C, N, O, F, P, S, Cl, Br and I. It rejects metals, radicals, query chemistry and virtual isotopic additions. Existing isotopic atoms are retained.

Every bond needs explicit aromaticity. A nonaromatic bond needs a declared integer order of 1, 2 or 3. An explicitly aromatic bond can retain fractional order 1.5 without an integer order, including a table consisting entirely of aromatic bonds. You do not need to invent alternating Kekule orders before adding H. Existing aromatic assignments and the chosen chemical state are preserved.

Native-to-RDKit conversion retains a declared prohibition of implicit H after sanitation, including aromatic NH. Forbidden inferred H still cause failure. The same prepared-component boundary covers the explicitly chosen 1QKU peptide fragment in the preparation recipe; it does not supply polymer boundary chemistry or repair missing heavy atoms.

The inventory tool separates existing indexed H atoms from virtual H annotations. Unknown counts are unresolved, not zero. Sanitation that changes the inventory or declared chemistry causes an error. Existing H are not removed or repositioned, including an excessive or conflicting input that must instead be rejected. Absolute R/S and E/Z assignments, when present, must agree with the supplied pose; unsupported descriptors need explicit normalization before this operation.

Requesting placement and reviewing the result#

After preparing molsys, call:

result = msm.build.add_missing_hydrogens(
    molsys, mode='fixed_chemical_state', pH=None, engine='RDKit',
    chemical_state='reference', structure_indices=[0], return_report=True,
)
molsys_with_h = result['molecular_system']
report = result['report']

The molsysmt.hydrogen_addition@1 report retains the state/structure indices, state ID, chemical readiness, per-parent inventory, generated IDs, original producer/engine versions, method, parameters, references and portable scientific attribution. atom_correspondence maps original to output indices, shape (n_original_atoms, 2); parent_hydrogen_pairs maps existing parents to appended H, shape (n_added_H, 2). The latter is empty with shape (0, 2) when no addition is needed. The report’s units are explicit; native coordinate storage is in nm. Stored virtual H counts become zero after their atoms have been materialized. Original atom indices, IDs, isotope assignments, memberships, coordinate values, box, structure/state links and existing bond chemistry are retained on a detached copy.

This method is local_hydrogen_placement, provided by RDKit AddHs(addCoords=True). Only the new coordinates come from that engine. Existing positions are copied from the source rather than roundtripped through its angstrom protocol.

Understanding validity and limits#

This operation does not select pH, protomers, tautomers or heavy-atom conformers, and does not optimize an energy or receptor-facing OH/NH orientation. H positions are generated local geometry, not experimentally observed coordinates. Separate refinement can be needed before a directional pharmacophore or docking analysis.

Supply an already compact periodic ligand. A bond requiring another periodic image is rejected; reconstruct explicitly with msm.pbc.wrap_to_mic() first. The source coordinates and box remain unchanged by H addition.

With attribute_policy='intersection', unsupported atom-aligned attributes (velocities, B-factors and occupancies), system observables and force-field parameters are dropped with a diagnostic and named in the report. New values are not fabricated. attribute_policy='strict' rejects such inputs. Sparse alternate locations for old atoms remain attached to their unchanged indices. Named interactions retain provenance but become unevaluated on an expanded output: adding an H can change interactions even if the old coordinates did not move. Recalculate explicitly. A zero-addition copy retains its analyses and attributes. Failure leaves the source unchanged.

Saving and retaining provenance#

Save molsys_with_h using the public H5MSM 0.5 conversion. It preserves the expanded chemical/coordinate domains, invalidated named analyses and historical preparation reports, even with return_report=False. Inspect independent copies through molsys_with_h.chemical_states.get_preparation_history(). Successful no-addition calls also append evidence; failures leave source history unchanged.

The final two records describe the general terminal-attachment operation and the RDKit H-generation operation that used it. Attachment records supplied coordinates; generation records local generated geometry, its original MolSysMT/RDKit versions, inventory and readiness audit. Nullable/mixed audit values are scalar lists; indices, masks and homogeneous columns retain explicit NumPy types. Both reports declare the original source atom/structure counts and evaluated structure_indices as an int64 array. They contain no coordinate snapshot.

Each envelope declares index_scope='operation': extraction, reordering or later edits do not remap historical indices or certify the current system. A state’s provenance index is not a pointer into this history. Retain intervening maps and caller template/input provenance separately. To carry an isolated component’s generation evidence into a larger destination, explicitly use destination.chemical_states.append_preparation_history( molsys_with_h.chemical_states.get_preparation_history()[-1:]). This archives evidence; it does not apply chemistry or align axes. See the reinsertion recipe.

The report contains NumPy arrays. Use the typed ChemicalStatesDict or H5MSM route, rather than assuming the dictionary is directly JSON serializable.

The existing default mode='pH', pH=7.4, engine='OpenMM' retains the legacy residue-oriented behavior. RDKit and report/state options require the explicit fixed-state mode. Neither mode predicts environment-dependent pKa.