MoleculeKit

repository·main·Indexed 42 days ago

https://github.com/acellera/moleculekit

A high-performance molecule reading, writing, and manipulation library designed for molecular modeling and analysis in cloud and web environments. It features a VMD-inspired atom-selection language, support for WebAssembly (Pyodide), and tools for voxelizing small molecules via SmallMolStack and generating tautomers with SmallMol.

Tokens
81.1K
Snippets
303
Records
371
Agent score
53%

What's inside MoleculeKit

  1. MoleculeKit How-to Guide Index

    main

    MoleculeKit provides a collection of task-oriented recipes for molecular modeling and analysis. The guides are organized into the following functional areas:

    • I/O: Reading and writing structures and trajectories.
    • Selection and manipulation: Atom selection, filtering, merging molecules, and managing properties.
    • Bonds, segments, topology: Guessing bonds, assigning segments/chains, and wrapping trajectories.
    • Geometry and analysis: Alignment, distance/contact computation, dihedrals, RMSD, and RMSF.
    • Projections and interactions: Computing projections and protein-ligand interactions.
    • Conversion to other toolkits: Converting data to RDKit and OpenFF formats.
    • Visualization: Viewing molecules, choosing representations, and rendering images.
    • External resources: Fetching data from RCSB and OPM.
  2. Explore the MoleculeKit API reference

    main

    The MoleculeKit API is organized into several functional domains to support molecular modeling, analysis, and visualization workflows. Key modules include:

    • Core: Fundamental objects like molecule, residues, readers, writers, and periodictable.
    • Atom Selection: Tools for selecting atoms using atomselect and related analysis utilities.
    • Geometry and Analysis: Functions for align, distance, dihedral, bondguesser, and unitcell operations.
    • Projections: A wide range of metric projections including metricrmsd, metricgyration, metricsecondarystructure, and more.
    • Interactions: Detection of molecular interactions such as hbonds, pipi, cationpi, and sigmahole.
    • File Formats: Support for various formats including pdbx (CIF), netcdf, dcd, trr, and xtc.
    • Tools: High-level utilities for preparation, mutation, modelling, docking, and obabel_tools.
    • Visualization: Integration with viewer, molstar, vmdviewer, and pymolgraphics.
    • Small Molecules: Specialized modules for small molecule handling via smallmol.
  3. Understand the KDTree implementation in MoleculeKit

    main

    MoleculeKit includes its own KDTree implementation derived from SciPy to avoid a direct dependency on the SciPy library.

    While it maintains high compatibility with SciPy's spatial algorithms, there are specific differences in the cKDTree.sparse_distance_matrix method:

    • Supported Output Types: Only 'dict' and 'ndarray' are supported.
    • Unsupported Types: The coo_array, coo_matrix, dok_array, and dok_matrix output types are not supported. If your workflow requires these specific sparse matrix formats, you must build them externally using the supported outputs.
  4. Control molecule drawing with Representations

    main

    Molecule drawing is controlled via mol.reps, which is a moleculekit.representations.Representations object. You can use mol.reps.add() to layer new visual styles onto your molecule or mol.reps.render() to output the result.

    Important Concept: Layering and Defaults Representations layer on top of each other. However, by default, adding a new representation replaces the entire automatic scene. To build upon the default automatic scene (which includes things like cartoon for proteins and waters), you must first call mol.reps.addDefaults(). This converts the automatic scene into editable entries.

    If you do not call addDefaults(), calling add() will wipe out all automatic representations (like the protein backbone) and only show what you explicitly added.

    from moleculekit.molecule import Molecule
    
    mol = Molecule("3PTB")
    
    # Start with the automatic scene as editable entries
    mol.reps.addDefaults()
    
    # Remove specific automatic entries (e.g., index 1 might be waters)
    mol.reps.remove(1)
    
    # Add custom layers
    mol.reps.add("protein and within 5 of resname BEN", "Surf", "Hydrophobicity", opacity=0.3)
    mol.reps.add("resname BEN", "ball-and-stick", "element-symbol", size=0.4)
    
    mol.render("pocket.png", center="resname BEN", zoom=0.3)
  5. Understand the Molecule data model

    main

    A moleculekit.molecule.Molecule object integrates two conceptual layers:

    1. Topology layer: Describes the atoms and their connections (e.g., name, resname, resid, chain, bonds). The number of atoms (numAtoms) is defined here.
    2. Trajectory layer: Carries per-frame data (e.g., coords, box, boxangles, fileloc, step, time). The number of frames (numFrames) is defined here.

    Standard workflow involves loading a topology file first, then appending trajectory files using mol.read(filename, append=True). Note that trajectories must have the same atom count as the topology, or an error will be raised.

  6. Use boolean masks for faster property access

    main

    If you already have a boolean mask (a numpy.ndarray), you can use direct array indexing instead of get and set. This is more concise because it avoids the string-based field lookup.

    • Reading: mol.field[mask] is equivalent to mol.get("field", sel=mask).
    • Writing: mol.field[mask] = value is equivalent to mol.set("field", value, sel=mask).

    Note on Coordinates:

    • mol.get("coords", sel=mask) returns the coordinates for the current active frame only.
    • mol.coords[mask] returns the full trajectory slice with shape (n_selection, 3, numFrames).
    # Assume ala_mask is a boolean mask for Alanine residues
    
    # Reading (Direct indexing is equivalent to .get)
    ala_names = mol.name[ala_mask]
    
    # Writing (Direct assignment is equivalent to .set)
    mol.charge[ala_mask] = 0.0
  7. Choose the correct file format for writing structures

    main

    Selecting a format depends on whether you need to preserve topology or just coordinates:

    • mmCIF (.cif, .bcif): Recommended default. Preserves full topology including bonds, bond orders, formal charges, segid, occupancy, and B-factors.
    • PDB (.pdb): Best for quick interchange with legacy tools. Note: Does NOT store bond orders. segid may be lost in many third-party parsers.
    • JSON (.json): Provides a lossless round-trip of the in-memory Molecule representation. Ideal for application state/development.
    • Trajectory Formats (.xtc, .dcd, .trr, .netcdf, .binpos): These write coordinates only. You must keep a topology file (like PSF or PDB) alongside them.
    • Small Molecules: Use .mol2 or .sdf to ensure bond orders are preserved.
    • Topology Only: Use .psf (CHARMM/NAMD).
  8. How the Molstar backend behaves

    main

    The molstar viewer is a built-in, browser-based backend that requires no external installation.

    Key Behaviors

    • Server Lifecycle: The first call to mol.view(viewer="molstar") starts an HTTP server on the first free port from 8765 onward and opens a browser tab. Subsequent calls to .view() for different molecules will push content to the same existing tab.
    • Automatic Updates: MoleculeKit uses a daemon thread to poll registered molecules every ~0.5s. If you mutate a molecule in Python, the changes propagate to the viewer automatically via Server-Sent Events; you do not need to call .view() again.
    • Mutation Handling:
      • Topology changes (e.g., mutation, filtering, bond edits) trigger a full rebuild of the scene.
      • Coordinate-only changes update the view in place without resetting the camera.
    • Representations: Representations set via mol.reps apply to the Molstar viewer. If a representation's selection matches no atoms, it is dropped with a warning. If all representations match nothing, view() will raise an error.
    • Multi-molecule Support: Multiple registered molecules appear in a slot sidebar (bottom-right) where they can be shown/hidden or removed. Removing a slot in the viewer does not affect the Python Molecule object.

    Troubleshooting Molstar

    • Headless Environments: Since Molstar opens a tab in the default web browser, it cannot automatically open a tab in headless environments.
    • Process Restarts: Restarting your Python process invalidates the existing browser tab. If you see a "Server restarted — refresh page" banner, reload the tab to reconnect.
    • Large Trajectories: For very large trajectories (>~100 MB of coordinates), consider striding the frames manually before calling .view() to avoid performance issues.
  9. Understanding Periodic Boundary Conditions (PBC) and Wrapping

    main

    In MD simulations, Periodic Boundary Conditions (PBC) are used to approximate infinite systems by treating the simulation box as if it were surrounded by identical copies (periodic images).

    Why molecules 'leave' the box

    MD engines often store absolute coordinates rather than resetting them to the central box every time an atom crosses a boundary. This prevents precision issues and allows for accurate measurement of total atomic drift. An atom 'leaving' the box is not a sign of a broken simulation; it has simply moved into a neighboring periodic image.

    What Wrapping does

    Wrapping is a visual and organizational tool. It does not change the underlying physics or the simulation results; it simply re-images the coordinates so that molecules appear centered and whole for easier visual inspection and distance calculations. It prevents molecules from appearing 'split' across the edges of the simulation box.

  10. How non-standard residue specs work

    main

    MoleculeKit uses four specific subclasses of spec objects to describe different types of non-standard residues. Understanding these helps you filter or manipulate the preparation process.

    1. ChainResidueSpec

    Used for residues inside a polypeptide chain that need special parameterization. This includes:

    • Non-canonical amino acids (NCAAs) embedded in a chain.
    • Canonical amino acids with sidechains covalently bonded to something else (e.g., N-glycosylated Asn, Glu–Lys isopeptide bonds).

    Attributes:

    • resname: Original residue name.
    • residue: UniqueResidueID (chain, resid, segid, insertion).
    • new_resname: The name to rename to before parameterization (None if no rename is needed).
    • anchor_atom: The atom involved in the non-peptide bond (None for plain NCAAs).
    • is_n_term / is_c_term: Boolean indicating if it is at the N- or C-terminus.

    2. CovalentLigandSpec

    Used for a residue that is not part of a chain but has exactly one covalent bond to the rest of the structure (e.g., a single-anchor covalent ligand).

    Attributes: resname, residue.

    3. LigandSpec

    Used for non-chain-resident residues with no covalent bonds to any other residue (free non-covalent ligands).

    Attributes: resname, residue.

    4. ScaffoldSpec

    Used for non-chain-resident residues with two or more covalent bonds to the polypeptide chain (e.g., bicyclic peptide scaffolds or multi-anchor covalent inhibitors).

    Attributes: resname, residue.

  11. Understand the internal layout of the coordinate array

    main

    MoleculeKit represents molecular dynamics trajectories as a three-dimensional NumPy array stored in mol.coords. This array has a shape of (numAtoms, 3, numFrames) using float32 precision.

    Array Axes:

    • Axis 0 (Atoms): Indexed 0 to mol.numAtoms - 1.
    • Axis 1 (Cartesian): The [x, y, z] coordinates in Ångström.
    • Axis 2 (Frames): Indexed 0 to mol.numFrames - 1.

    Key Properties:

    • mol.frame: An integer tracking the current frame index. This index is used by the viewer and for distance-based selections like within or exwithin. You can update it using mol.frame = i.
    • For single-frame structures, mol.numFrames == 1.
    # Common slicing patterns for mol.coords (shape: numAtoms, 3, numFrames)
    
    # All atoms, frame 0: shape (numAtoms, 3)
    frame_0 = mol.coords[:, :, 0] 
    
    # All frames for atom 42: shape (3, numFrames)
    atom_42 = mol.coords[42, :, :]
    
    # x-coordinates of all atoms, all frames: (numAtoms, numFrames)
    x_coords = mol.coords[:, 0, :]
    
    # The last frame
    last_frame = mol.coords[:, :, -1]