MDAnalysis - Python Toolkit for Molecular Dynamics Analysis

pip install MDAnalysis

import MDAnalysis as mda

# Core classes available at top level
from MDAnalysis import Universe, AtomGroup, ResidueGroup, SegmentGroup, Writer
from MDAnalysis import Merge  # For merging universes

# Exception classes
from MDAnalysis.exceptions import (
    SelectionError, NoDataError, ApplicationError,
    SelectionWarning, MissingDataWarning
)

# Units handling
import MDAnalysis.units

# Analysis modules (imported as needed)
import MDAnalysis.analysis.align
import MDAnalysis.analysis.rms
import MDAnalysis.analysis.distances

def __init__(self, topology, *coordinates, **kwargs):
    """
    Create a Universe from topology and coordinate files.
    
    Parameters
    ----------
    topology : str or file-like
        Topology file (PSF, PDB, GRO, TPR, etc.) containing atomic structure
    *coordinates : str or file-like, optional
        Trajectory file(s) (DCD, XTC, TRR, etc.) with coordinate time series
    **kwargs : dict
        Additional keyword arguments passed to readers
        
    Returns
    -------
    Universe
        Universe object ready for analysis
    """

# Basic usage examples
u = mda.Universe("topology.psf", "trajectory.dcd")
u = mda.Universe("structure.pdb")  # Single frame from PDB
u = mda.Universe("topol.tpr", "traj.xtc", "traj2.xtc")  # Multiple trajectories

# Basic selections
protein = u.select_atoms("protein")
ca_atoms = u.select_atoms("name CA")
backbone = u.select_atoms("backbone")

# Residue-based selections  
active_site = u.select_atoms("resid 23 45 67")
binding_pocket = u.select_atoms("resname HIS PHE TRP")

# Geometric selections
near_protein = u.select_atoms("around 5.0 protein")
water_shell = u.select_atoms("resname SOL and around 8.0 protein")

# Boolean operations
flexible_region = u.select_atoms("protein and not (name CA or name C or name N)")

# Iterate through trajectory frames
for ts in u.trajectory:
    print(f"Frame {ts.frame}, Time: {ts.time} ps")
    
    # Perform per-frame analysis
    com = protein.center_of_mass()
    rgyr = protein.radius_of_gyration()
    
# Random access to specific frames
u.trajectory[0]    # First frame
u.trajectory[-1]   # Last frame
u.trajectory[100:200:10]  # Slice with step

# Universe operations
u = mda.Universe("topology.pdb", "trajectory.xtc") 
merged = mda.Merge(u1.atoms, u2.atoms)  # Combine systems

# Selection and grouping
protein = u.select_atoms("protein")
ca_atoms = protein.select_atoms("name CA")
binding_site = u.select_atoms("resid 23-45 and around 5.0 resname LIG")

# Basic analysis
center = protein.center_of_mass()
radius = protein.radius_of_gyration() 
coords = protein.positions  # Current frame coordinates

# Reading various formats
u1 = mda.Universe("system.psf", "trajectory.dcd")  # CHARMM
u2 = mda.Universe("topol.tpr", "traj.xtc")        # GROMACS  
u3 = mda.Universe("structure.prmtop", "mdcrd")     # AMBER

# Writing coordinate data
protein.write("protein_only.pdb")
u.atoms.write("all_atoms.gro")

# Trajectory writing
with mda.Writer("output.dcd", u.atoms.n_atoms) as W:
    for ts in u.trajectory:
        W.write(u.atoms)

from MDAnalysis.analysis import align, rms, distances

# RMSD analysis
R = rms.RMSD(u, u, select="backbone")
R.run()
print(R.results.rmsd)  # Time series of RMSD values

# Distance analysis
dist_array = distances.distance_array(group1.positions, group2.positions)
contacts = distances.contact_matrix(protein.positions, cutoff=8.0)

# Structural alignment
alignto_result = align.alignto(mobile, reference, select="name CA")

# Topology queries
print(u.atoms.bonds)     # Bond connectivity
print(u.atoms.angles)    # Angle definitions
print(u.atoms.dihedrals) # Dihedral angles

# Adding topology information
u.add_bonds([(0, 1), (1, 2), (2, 3)])  # Add bonds by atom indices
u.guess_bonds()  # Automatically guess bonds from distances

# Topology attributes
u.add_TopologyAttr('masses', [12.0] * u.atoms.n_atoms)
masses = u.atoms.masses

# Basic transformations
protein.translate([10.0, 0.0, 0.0])  # Move along x-axis
protein.rotate([[1,0,0], [0,0,1], [0,1,0]])  # Rotation matrix
protein.transform(transformation_matrix)  # General transformation

# Alignment operations  
align.alignto(mobile, reference, select="name CA")
transformations.fit_rot_trans(mobile, reference)

# Periodic boundary handling
protein.wrap(compound="residues")      # Wrap by residue
protein.unwrap(compound="fragments")   # Unwrap fragments

# Property-based selections
u.select_atoms("name CA CB CG")
u.select_atoms("resname ALA VAL LEU ILE")  
u.select_atoms("resid 1-100 and name CA")

# Geometric selections
u.select_atoms("around 5.0 resname LIG")          # Within distance
u.select_atoms("spherical_layer 2.0 4.0 protein") # Spherical shell
u.select_atoms("sphlayer 2.0 4.0 protein")        # Alias for above

# Complex boolean selections  
u.select_atoms("protein and not backbone")
u.select_atoms("(resname ALA or resname GLY) and name CA")

import MDAnalysis.units as units

# Unit conversions
units.convert(10.0, 'Angstrom', 'nm')  # Length conversion
units.convert(300.0, 'K', 'kJ/mol')    # Temperature to energy

# Base units in MDAnalysis
print(units.MDANALYSIS_BASE_UNITS)
# {'length': 'Angstrom', 'time': 'ps', 'energy': 'kJ/mol', ...}

# Distance calculations
from MDAnalysis.lib.distances import distance_array, self_distance_array
distances = distance_array(coords1, coords2, box=u.dimensions)

from MDAnalysis.auxiliary import auxreader, get_auxreader_for
from MDAnalysis.auxiliary import XVGReader, EDRReader

# Load auxiliary data
aux = auxreader('pullforce.xvg')
u.trajectory.add_auxiliary('force', aux)

# Access during iteration
for ts in u.trajectory:
    force_value = ts.aux.force

from MDAnalysis.converters import RDKit, ParmEd, OpenMM

# Convert to RDKit molecule
rdkit_mol = protein.convert_to("RDKIT")

# Convert to ParmEd structure  
parmed_struct = u.atoms.convert_to("PARMED")

# Convert to OpenMM topology
topology, positions = u.atoms.convert_to("OPENMM")