molecular-dynamics

Molekulardynamik

Übersicht

Molekulardynamik (MD)-Simulationen modellieren rechnerisch die zeitliche Entwicklung molekularer Systeme durch Integration der Newtonschen Bewegungsgleichungen. Dieser Skill deckt zwei komplementäre Werkzeuge ab:

• OpenMM (https://openmm.org/): Hochleistungs-MD-Simulationsengine mit GPU-Unterstützung, Python-API und flexibler Kraftfeldunterstützung
• MDAnalysis (https://mdanalysis.org/): Python-Bibliothek zum Lesen, Schreiben und Analysieren von MD-Trajektorien aus allen gängigen Simulationspaketen

Installation:

conda install -c conda-forge openmm mdanalysis nglview
# oder
uv pip install openmm mdanalysis

Wann dieser Skill verwendet werden sollte

Verwenden Sie Molekulardynamik, wenn:

• Proteinstabilitätsanalyse: Wie wirkt sich eine Mutation auf die Proteindynamik aus?
• Wirkstoffbindungssimulationen: Charakterisieren Sie den Bindungsmodus und die Verweilzeit eines Liganden
• Konformationsabtastung: Erkunden Sie die Proteinflexibilität und Konformationsänderungen
• Protein-Protein-Wechselwirkung: Modellieren Sie Grenzflächendynamik und Bindungsenergetik
• RMSD/RMSF-Analyse: Quantifizieren Sie strukturelle Fluktuationen aus einer Referenzstruktur
• Freie-Energie-Schätzung: Berechnen Sie die freie Bindungsenergie oder die freie Konformationsenergie
• Membransimulationen: Modellieren Sie Proteine in Lipiddoppelschichten
• Intrinsisch ungeordnete Proteine: Untersuchen Sie konformationelle Ensembles von IDPs

Kern-Workflow: OpenMM-Simulation

1. Systemvorbereitung

from openmm.app import *
from openmm import *
from openmm.unit import *
import sys

def prepare_system_from_pdb(pdb_datei, kraftfeld_name="amber14-all.xml",
                              wasser_modell="amber14/tip3pfb.xml"):
    """
    Bereitet ein OpenMM-System aus einer PDB-Datei vor.

    Args:
        pdb_datei: Pfad zur bereinigten PDB-Datei (verwenden Sie PDBFixer für rohe PDB-Dateien)
        kraftfeld_name: Kraftfeld-XML-Datei
        wasser_modell: Wassermodell-XML-Datei

    Returns:
        pdb, kraftfeld, system, topologie
    """
    # PDB laden
    pdb = PDBFile(pdb_datei)

    # Kraftfeld laden
    kraftfeld = ForceField(kraftfeld_name, wasser_modell)

    # Wasserstoffatome hinzufügen und solvatisieren
    modellierer = Modeller(pdb.topology, pdb.positions)
    modellierer.addHydrogens(kraftfeld)

    # Lösungsmittelbox hinzufügen (10 Å Abstand, 150 mM NaCl)
    modellierer.addSolvent(
        kraftfeld,
        model='tip3p',
        padding=10*angstroms,
        ionicStrength=0.15*molar
    )

    print(f"System: {modellierer.topology.getNumAtoms()} Atome, "
          f"{modellierer.topology.getNumResidues()} Reste")

    # System erstellen
    system = kraftfeld.createSystem(
        modellierer.topology,
        nonbondedMethod=PME,         # Particle Mesh Ewald für langreichweitige Elektrostatik
        nonbondedCutoff=1.0*nanometer,
        constraints=HBonds,           # Wasserstoffbrückenbindungen einschränken (ermöglicht 2-fs-Schrittweite)
        rigidWater=True,
        ewaldErrorTolerance=0.0005
    )

    return modellierer, system

2. Energieminimierung

from openmm.app import *
from openmm import *
from openmm.unit import *

def minimize_energy(modellierer, system, ausgabe_pdb="minimiert.pdb",
                     max_iterationen=1000, toleranz=10.0):
    """
    Energieminimierung des Systems, um sterische Hinderungen zu beseitigen.

    Args:
        modellierer: Modeller-Objekt mit Topologie und Positionen
        system: OpenMM-System
        ausgabe_pdb: Pfad zum Speichern der minimierten Struktur
        max_iterationen: Maximale Anzahl von Minimierungsschritten
        toleranz: Konvergenzkriterium in kJ/mol/nm

    Returns:
        Simulationsobjekt mit minimierten Positionen
    """
    # Integrator einrichten (spielt für die Minimierung keine Rolle)
    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)

    # Simulation erstellen
    # GPU verwenden, falls verfügbar (CUDA oder OpenCL), sonst auf CPU zurückgreifen
    try:
        platform = Platform.getPlatformByName('CUDA')
        eigenschaften = {'DeviceIndex': '0', 'Precision': 'mixed'}
    except Exception:
        try:
            platform = Platform.getPlatformByName('OpenCL')
            eigenschaften = {}
        except Exception:
            platform = Platform.getPlatformByName('CPU')
            eigenschaften = {}

    simulation = Simulation(
        modellierer.topology, system, integrator,
        platform, eigenschaften
    )
    simulation.context.setPositions(modellierer.positions)

    # Anfangsenergie überprüfen
    zustand = simulation.context.getState(getEnergy=True)
    print(f"Anfangsenergie: {zustand.getPotentialEnergy()}")

    # Minimieren
    simulation.minimizeEnergy(
        tolerance=toleranz*kilojoules_per_mole/nanometer,
        maxIterations=max_iterationen
    )

    zustand = simulation.context.getState(getEnergy=True, getPositions=True)
    print(f"Minimierte Energie: {zustand.getPotentialEnergy()}")

    # Minimierte Struktur speichern
    with open(ausgabe_pdb, 'w') as f:
        PDBFile.writeFile(simulation.topology, zustand.getPositions(), f)

    return simulation

3. NVT-Äquilibrierung

from openmm.app import *
from openmm import *
from openmm.unit import *

def run_nvt_aequilibrierung(simulation, n_schritte=50000, temperatur=300,
                            berichtsintervall=1000, ausgabe_praefix="nvt"):
    """
    NVT-Äquilibrierung: konstantes N, V, T.
    Äquilibriert Geschwindigkeiten auf Zieltemperatur.

    Args:
        simulation: OpenMM-Simulation (nach Minimierung)
        n_schritte: Anzahl der MD-Schritte (50000 × 2fs = 100 ps)
        temperatur: Temperatur in Kelvin
        berichtsintervall: Schritte zwischen Datenberichten
        ausgabe_praefix: Dateipräfix für Trajektorie und Log
    """
    # Positionsbeschränkungen für Rückgrat während NVT hinzufügen
    # (Optional: schwere Atome einschränken)

    # Temperatur einstellen
    simulation.context.setVelocitiesToTemperature(temperatur*kelvin)

    # Berichterstatter hinzufügen
    simulation.reporters = []

    # Logdatei
    simulation.reporters.append(
        StateDataReporter(
            f"{ausgabe_praefix}_log.txt",
            berichtsintervall,
            step=True,
            potentialEnergy=True,
            kineticEnergy=True,
            temperature=True,
            volume=True,
            speed=True
        )
    )

    # DCD-Trajektorie (kompaktes Binärformat)
    simulation.reporters.append(
        DCDReporter(f"{ausgabe_praefix}_traj.dcd", berichtsintervall)
    )

    print(f"Führe NVT-Äquilibrierung durch: {n_schritte} Schritte ({n_schritte*2/1000:.1f} ps)")
    simulation.step(n_schritte)
    print("NVT-Äquilibrierung abgeschlossen")

    return simulation

4. NPT-Äquilibrierung und Produktion

def run_npt_produktion(simulation, n_schritte=500000, temperatur=300, druck=1.0,
                        berichtsintervall=5000, ausgabe_praefix="npt"):
    """
    NPT-Produktionslauf: konstantes N, P, T.

    Args:
        n_schritte: Produktionsschritte (500000 × 2fs = 1 ns)
        temperatur: Temperatur in Kelvin
        druck: Druck in bar
        berichtsintervall: Schritte zwischen Berichten
    """
    # Monte-Carlo-Barostat zur Druckkontrolle hinzufügen
    system = simulation.context.getSystem()
    system.addForce(MonteCarloBarostat(druck*bar, temperatur*kelvin, 25))
    simulation.context.reinitialize(preserveState=True)

    # Berichterstatter aktualisieren
    simulation.reporters = []
    simulation.reporters.append(
        StateDataReporter(
            f"{ausgabe_praefix}_log.txt",
            berichtsintervall,
            step=True,
            potentialEnergy=True,
            temperature=True,
            density=True,
            speed=True
        )
    )
    simulation.reporters.append(
        DCDReporter(f"{ausgabe_praefix}_traj.dcd", berichtsintervall)
    )

    # Prüfpunkte speichern
    simulation.reporters.append(
        CheckpointReporter(f"{ausgabe_praefix}_checkpoint.chk", 50000)
    )

    print(f"Führe NPT-Produktion durch: {n_schritte} Schritte ({n_schritte*2/1000000:.2f} ns)")
    simulation.step(n_schritte)
    print("Produktions-MD abgeschlossen")
    return simulation

Trajektorienanalyse mit MDAnalysis

1. Trajektorie laden

import MDAnalysis as mda
from MDAnalysis.analysis import rms, align, contacts
import numpy as np
import matplotlib.pyplot as plt

def load_trajectory(topologie_datei, trajektorien_datei):
    """
    Lädt eine MD-Trajektorie mit MDAnalysis.

    Args:
        topologie_datei: PDB-, PSF- oder andere Topologiedatei
        trajektorien_datei: DCD-, XTC-, TRR- oder andere Trajektorie
    """
    u = mda.Universe(topologie_datei, trajektorien_datei)
    print(f"Universum: {u.atoms.n_atoms} Atome, {u.trajectory.n_frames} Frames")
    print(f"Zeitspanne: 0 bis {u.trajectory.totaltime:.0f} ps")
    return u

2. RMSD-Analyse

def compute_rmsd(u, selection="backbone", referenz_frame=0):
    """
    Berechnet den RMSD ausgewählter Atome relativ zum Referenzframe.

    Args:
        u: MDAnalysis-Universum
        selection: Atomauswahlzeichenfolge (MDAnalysis-Syntax)
        referenz_frame: Frame-Index für die Referenzstruktur

    Returns:
        numpy-Array mit (Zeit, RMSD)-Werten
    """
    # Trajektorie ausrichten, um RMSD zu minimieren
    aligner = align.AlignTraj(u, u, select=selection, in_memory=True)
    aligner.run()

    # RMSD berechnen
    R = rms.RMSD(u, select=selection, ref_frame=referenz_frame)
    R.run()

    rmsd_daten = R.results.rmsd  # Spalten: Frame, Zeit, RMSD
    return rmsd_daten

def plot_rmsd(rmsd_daten, titel="RMSD über die Zeit", ausgabedatei="rmsd.png"):
    """Zeichnet RMSD über die Simulationszeit."""
    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(rmsd_daten[:, 1] / 1000, rmsd_daten[:, 2], 'b-', linewidth=0.5)
    ax.set_xlabel("Zeit (ns)")
    ax.set_ylabel("RMSD (Å)")
    ax.set_title(titel)
    ax.axhline(rmsd_daten[:, 2].mean(), color='r', linestyle='--',
               label=f'Mittelwert: {rmsd_daten[:, 2].mean():.2f} Å')
    ax.legend()
    plt.tight_layout()
    plt.savefig(ausgabedatei, dpi=150)
    return fig

3. RMSF-Analyse (Flexibilität pro Rest)

def compute_rmsf(u, selection="backbone", start_frame=0):
    """
    Berechnet die RMSF pro Rest (Flexibilität).

    Returns:
        residue_ids, rmsf_werte Arrays
    """
    # Atome auswählen
    atome = u.select_atoms(selection)

    # RMSF berechnen
    R = rms.RMSF(atome)
    R.run(start=start_frame)

    # Nach Rest mitteln
    residue_ids = []
    rmsf_pro_rest = []
    for res in u.select_atoms(selection).residues:
        res_atome = res.atoms.intersection(atome)
        if len(res_atome) > 0:
            residue_ids.append(res.resid)
            rmsf_pro_rest.append(R.results.rmsf[res_atome.indices].mean())

    return np.array(residue_ids), np.array(rmsf_pro_rest)

4. Protein-Ligand-Kontakte

def analyze_contacts(u, protein_sel="protein", ligand_sel="resname LIG",
                      radius=4.5, start_frame=0):
    """
    Verfolgt Protein-Ligand-Kontakte über die Trajektorie.

    Args:
        radius: Kontaktdistanz-Cutoff in Angstrom
    """
    protein = u.select_atoms(protein_sel)
    ligand = u.select_atoms(ligand_sel)

    kontakt_frames = []
    for ts in u.trajectory[start_frame:]:
        # Proteinatome innerhalb des Radius um den Liganden finden
        distanzen = contacts.contact_matrix(
            protein.positions, ligand.positions, radius
        )
        kontakt_reste = set()
        for i in range(distanzen.shape[0]):
            if distanzen[i].any():
                kontakt_reste.add(protein.atoms[i].resid)
        kontakt_frames.append(kontakt_reste)

    return kontakt_frames

Kraftfeld-Auswahlhilfe

System	Empfohlenes Kraftfeld	Wassermodell
Standardproteine	AMBER14 (amber14-all.xml)	TIP3P-FB
Proteine + kleine Moleküle	AMBER14 + GAFF2	TIP3P-FB
Membranproteine	CHARMM36m	TIP3P
Nukleinsäuren	AMBER99-bsc1 oder AMBER14	TIP3P
Ungeordnete Proteine	ff19SB oder CHARMM36m	TIP3P

Werkzeuge zur Systemvorbereitung

PDBFixer (für rohe PDB-Dateien)

from pdbfixer import PDBFixer
from openmm.app import PDBFile

def fix_pdb(eingabe_pdb, ausgabe_pdb, ph=7.0):
    """Behebt häufige PDB-Probleme: fehlende Reste, Atome, fügt H hinzu, standardisiert."""
    fixer = PDBFixer(filename=eingabe_pdb)
    fixer.findMissingResidues()
    fixer.findNonstandardResidues()
    fixer.replaceNonstandardResidues()
    fixer.removeHeterogens(True)    # Wasser/Liganden entfernen
    fixer.findMissingAtoms()
    fixer.addMissingAtoms()
    fixer.addMissingHydrogens(ph)

    with open(ausgabe_pdb, 'w') as f:
        PDBFile.writeFile(fixer.topology, fixer.positions, f)

    return ausgabe_pdb

GAFF2 für kleine Moleküle (über OpenFF-Toolkit)

# Für die Ligandenparametrisierung das OpenFF-Toolkit oder ACPYPE verwenden
# uv pip install openff-toolkit
from openff.toolkit import Molecule, ForceField as OFFForceField
from openff.interchange import Interchange

def parameterize_ligand(smiles, ff_name="openff-2.0.0.offxml"):
    """Erzeugt GAFF2/OpenFF-Parameter für ein kleines Molekül."""
    mol = Molecule.from_smiles(smiles)
    mol.generate_conformers(n_conformers=1)

    off_ff = OFFForceField(ff_name)
    interchange = off_ff.create_interchange(mol.to_topology())
    return interchange

Best Practices

• Vor MD immer minimieren: Rohe PDB-Strukturen weisen sterische Hinderungen auf
• Vor der Produktion äquilibrieren: NVT (50–100 ps) → NPT (100–500 ps) → Produktion
• GPU verwenden: Simulationen sind auf GPU (CUDA/OpenCL) 10–100× schneller
• 2-fs-Schrittweite mit HBonds-Einschränkungen: Standard; 4 fs mit HMR (Wasserstoffmassen-Neuverteilung)
• Nur äquilibrierte Trajektorie analysieren: Erste 20–50 % als Äquilibrierung verwerfen
• Prüfpunkte speichern: MD-Läufe können fehlschlagen; Prüfpunkte ermöglichen Neustarts
• Periodische Randbedingungen: Für solvatisierte Systeme erforderlich
• PME für Elektrostatik: Genauer als Cutoff-Methoden für geladene Systeme

Zusätzliche Ressourcen

• OpenMM-Dokumentation: https://openmm.org/documentation.html
• MDAnalysis-Benutzerhandbuch: https://docs.mdanalysis.org/
• GROMACS (alternative MD-Engine): https://manual.gromacs.org/
• NAMD (Alternative): https://www.ks.uiuc.edu/Research/namd/
• CHARMM-GUI (Webbasierte Systemvorbereitung): https://www.charmm-gui.org/