molecular-dynamics

分子動力学

概要

分子動力学（MD）シミュレーションは、ニュートンの運動方程式を積分することで、分子系の時間発展を計算機上でモデル化します。このスキルでは、相補的な2つのツールを扱います。

• OpenMM（https://openmm.org/）: GPU サポート、Python API、柔軟な力場サポートを備えた高性能 MD シミュレーションエンジン
• MDAnalysis（https://mdanalysis.org/）: すべての主要なシミュレーションパッケージの MD トラジェクトリを読み書き・解析するための Python ライブラリ

レビュー済みバージョン: OpenMM 8.6.1 および MDAnalysis 2.10.0。参照プラットフォームでのスモークチェックは小規模な合成系を使用しており、GPU 性能や科学的な収束はそれらによって確立されるものではありません。OpenMM の最新 API ページは開発ビルドを示しているため、インストール済みの 8.6.1 API が実行可能なリファレンスとなります。

専用環境へのインストール:

uv venv --python 3.12 .venv-md
uv pip install --python .venv-md/bin/python "openmm==8.6.1" "MDAnalysis==2.10.0" matplotlib pandas
# Windows: use .venv-md/Scripts/python.exe instead

このスキルを使うべき場面

次のような場合に分子動力学を使用します。

• タンパク質安定性解析: 変異はタンパク質の動態にどのように影響するか？
• 薬物結合シミュレーション: リガンドの結合様式と滞留時間を特性評価する
• コンフォメーションサンプリング: タンパク質の柔軟性とコンフォメーション変化を探索する
• タンパク質間相互作用: 界面の動態と結合エネルギーをモデル化する
• RMSD/RMSF 解析: 参照構造からの構造揺らぎを定量化する
• 自由エネルギー推定: 結合自由エネルギーまたはコンフォメーション自由エネルギーを計算する
• 膜シミュレーション: 脂質二重層中のタンパク質をモデル化する
• 天然変性タンパク質: IDR のコンフォメーションアンサンブルを研究する

コアワークフロー: OpenMM シミュレーション

例には系固有の準備と検証が必要です。パラメータ化の前に、生体集合体、代替位置、欠損ループ、末端、プロトン化、ジスルフィド、リガンド、補因子を確認してください。MD トラジェクトリだけでは結合自由エネルギーや滞留時間は確立されません。

以下の関数は1つの Python モジュールを構成します。それらの import ブロックはまとめて実行してください。以前の結果を上書きしないよう、各段階で新しい出力プレフィックスを使用してください。

1. 系の準備

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

def prepare_system_from_pdb(pdb_file, forcefield_name="amber14-all.xml",
                              water_model="amber14/tip3pfb.xml", ph=7.0):
    """
    Prepare an OpenMM system from a PDB file.

    Args:
        pdb_file: Path to cleaned PDB file (use PDBFixer for raw PDB files)
        forcefield_name: Force field XML file
        water_model: Water model XML file

    Returns:
        modeller, system
    """
    # Load PDB
    pdb = PDBFile(pdb_file)

    # Load force field
    forcefield = ForceField(forcefield_name, water_model)

    # Add hydrogens and solvate
    modeller = Modeller(pdb.topology, pdb.positions)
    modeller.addHydrogens(forcefield, pH=ph)

    # TIP3P geometry also serves TIP3P-FB; XML supplies the parameters.
    # 10 Å padding; added salt excludes neutralizing counterions.
    modeller.addSolvent(
        forcefield,
        model='tip3p',
        padding=10*angstroms,
        ionicStrength=0.15*molar, neutralize=True
    )

    print(f"System: {modeller.topology.getNumAtoms()} atoms, "
          f"{modeller.topology.getNumResidues()} residues")

    # Create system
    system = forcefield.createSystem(
        modeller.topology,
        nonbondedMethod=PME,         # Particle Mesh Ewald for long-range electrostatics
        nonbondedCutoff=1.0*nanometer,
        constraints=HBonds,           # Constrain bonds involving H (not intermolecular H-bonds)
        rigidWater=True,
        ewaldErrorTolerance=0.0005
    )

    return modeller, system

2. エネルギー最小化

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

def minimize_energy(modeller, system, output_pdb="minimized.pdb",
                     max_iterations=1000, tolerance=10.0, platform_name=None):
    """
    Energy minimize the system to remove steric clashes.

    Args:
        modeller: Modeller object with topology and positions
        system: OpenMM System
        output_pdb: Path to save minimized structure
        max_iterations: Maximum minimization steps
        tolerance: Convergence criterion in kJ/mol/nm

    Returns:
        simulation object with minimized positions
    """
    # Keep the integrator at 2 fs for the subsequent NVT/NPT examples.
    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)

    # Automatic platform selection, or an explicit tested platform (e.g. CPU).
    # A registered GPU plugin does not prove a working device or driver.
    platform = Platform.getPlatformByName(platform_name) if platform_name else None
    simulation = Simulation(modeller.topology, system, integrator, platform)
    simulation.context.setPositions(modeller.positions)

    # Check initial energy
    state = simulation.context.getState(getEnergy=True)
    print(f"Initial energy: {state.getPotentialEnergy()}")

    # Minimize
    simulation.minimizeEnergy(
        tolerance=tolerance*kilojoules_per_mole/nanometer,
        maxIterations=max_iterations
    )

    state = simulation.context.getState(getEnergy=True, getPositions=True)
    print(f"Minimized energy: {state.getPotentialEnergy()}")

    # Save minimized structure
    with open(output_pdb, 'w') as f:
        PDBFile.writeFile(simulation.topology, state.getPositions(), f)

    return simulation

3. NVT 平衡化

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

def run_nvt_equilibration(simulation, n_steps=50000, temperature=300,
                            report_interval=1000, output_prefix="nvt"):
    """
    NVT equilibration: constant N, V, T.
    Equilibrate velocities to target temperature.

    Args:
        simulation: OpenMM Simulation (after minimization)
        n_steps: Number of MD steps (50000 × 2fs = 100 ps)
        temperature: Temperature in Kelvin
        report_interval: Steps between data reports
        output_prefix: File prefix for trajectory and log
    """
    # This example is unrestrained; add validated restraints before Context creation.
    if any("Barostat" in type(f).__name__ and f.getFrequency() > 0
           for f in simulation.system.getForces()):
        raise ValueError("NVT requires all barostats disabled or absent")

    # Set both the thermostat target and initial velocities.
    simulation.integrator.setTemperature(temperature*kelvin)
    simulation.context.setVelocitiesToTemperature(temperature*kelvin)

    # Add reporters
    simulation.reporters = []

    # Log file
    simulation.reporters.append(
        StateDataReporter(
            f"{output_prefix}_log.txt",
            report_interval,
            step=True,
            potentialEnergy=True,
            kineticEnergy=True,
            temperature=True,
            volume=True,
            speed=True
        )
    )

    # DCD trajectory (compact binary format)
    simulation.reporters.append(
        DCDReporter(f"{output_prefix}_traj.dcd", report_interval)
    )

    duration = (n_steps * simulation.integrator.getStepSize()).value_in_unit(picoseconds)
    print(f"Running NVT equilibration: {n_steps} steps ({duration:.1f} ps)")
    simulation.step(n_steps)
    print("NVT equilibration complete")

    return simulation

4. NPT 平衡化と本番計算

まずこの関数を専用の NPT 平衡化プレフィックスで呼び出してください。別途本番計算を呼び出す前に、密度、エネルギー、構造、および再現安定性を確認してください。デフォルトの時間は一例であり、平衡化や収束の基準ではありません。

def run_npt_production(simulation, n_steps=500000, temperature=300, pressure=1.0,
                        report_interval=5000, output_prefix="npt"):
    """
    NPT production run: constant N, P, T.

    Args:
        n_steps: Production steps (500000 × 2fs = 1 ns)
        temperature: Temperature in Kelvin
        pressure: Pressure in bar
        report_interval: Steps between reports
    """
    # Keep the Langevin thermostat and barostat at the same temperature.
    simulation.integrator.setTemperature(temperature*kelvin)
    system = simulation.system
    barostats = [f for f in system.getForces() if "Barostat" in type(f).__name__]
    if not system.usesPeriodicBoundaryConditions():
        raise ValueError("NPT requires a periodic system")
    if not barostats:
        system.addForce(MonteCarloBarostat(pressure*bar, temperature*kelvin, 25))
        simulation.context.reinitialize(preserveState=True)
    elif len(barostats) != 1 or not isinstance(barostats[0], MonteCarloBarostat):
        raise ValueError("This example supports one isotropic MonteCarloBarostat")
    else:
        barostats[0].setDefaultPressure(pressure*bar)
        barostats[0].setDefaultTemperature(temperature*kelvin)
        barostats[0].setFrequency(25)
    # Existing Context parameters must also be updated on repeated calls.
    simulation.context.setParameter(MonteCarloBarostat.Pressure(), pressure)
    simulation.context.setParameter(MonteCarloBarostat.Temperature(), temperature)

    # Update reporters
    simulation.reporters = []
    simulation.reporters.append(
        StateDataReporter(
            f"{output_prefix}_log.txt",
            report_interval,
            step=True,
            potentialEnergy=True,
            temperature=True,
            density=True,
            speed=True
        )
    )
    simulation.reporters.append(
        DCDReporter(f"{output_prefix}_traj.dcd", report_interval)
    )

    # Save checkpoints
    simulation.reporters.append(
        CheckpointReporter(f"{output_prefix}_checkpoint.chk", 50000)
    )

    duration = (n_steps * simulation.integrator.getStepSize()).value_in_unit(nanoseconds)
    print(f"Running NPT stage: {n_steps} steps ({duration:.3f} ns)")
    simulation.step(n_steps)
    simulation.saveCheckpoint(f"{output_prefix}_checkpoint.chk")
    simulation.saveState(f"{output_prefix}_state.xml")
    print("NPT stage complete")
    return simulation

MDAnalysis によるトラジェクトリ解析

1. トラジェクトリの読み込み

元の非溶媒和入力ではなく、DCD の原子数と順序に正確に対応する溶媒和済みトポロジー（たとえば上記の minimized.pdb）を使用してください。MDAnalysis はデフォルトで長さを Å、時間を ps に変換しますが、OpenMM の生の座標は nm を使用します。フレームのボックスベクトルと実際のタイムスタンプを保持してください。フレーム番号から時間を推測しないでください。

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

def load_trajectory(topology_file, trajectory_file):
    """
    Load an MD trajectory with MDAnalysis.

    Args:
        topology_file: PDB, PSF, or other topology file
        trajectory_file: DCD, XTC, TRR, or other trajectory
    """
    u = mda.Universe(topology_file, trajectory_file)
    print(f"Universe: {u.atoms.n_atoms} atoms, {u.trajectory.n_frames} frames")
    first, last = u.trajectory[0].time, u.trajectory[-1].time
    print(f"Time range: {first:g} to {last:g} ps")
    u.trajectory[0]
    return u

2. RMSD 解析

def compute_rmsd(u, selection="backbone", reference_frame=0):
    """
    Compute RMSD of selected atoms relative to reference frame.

    Args:
        u: MDAnalysis Universe
        selection: Atom selection string (MDAnalysis syntax)
        reference_frame: Frame index for reference structure

    Returns:
        numpy array with columns [frame index, time in ps, RMSD in Angstroms]
    """
    # RMSD performs its own fit; avoid rotating the stored trajectory here.
    # Make the selected molecule whole before analysis of periodic trajectories.
    R = rms.RMSD(u, select=selection, ref_frame=reference_frame)
    R.run()

    rmsd_data = R.results.rmsd  # columns: frame, time, RMSD
    return rmsd_data

def plot_rmsd(rmsd_data, title="RMSD over time", output_file="rmsd.png"):
    """Plot RMSD over simulation time."""
    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(rmsd_data[:, 1] / 1000, rmsd_data[:, 2], 'b-', linewidth=0.5)
    ax.set_xlabel("Time (ns)")
    ax.set_ylabel("RMSD (Å)")
    ax.set_title(title)
    ax.axhline(rmsd_data[:, 2].mean(), color='r', linestyle='--',
               label=f'Mean: {rmsd_data[:, 2].mean():.2f} Å')
    ax.legend()
    plt.tight_layout()
    plt.savefig(output_file, dpi=150)
    return fig

3. RMSF 解析（残基ごとの柔軟性）

RMSF の前に分子を完全なものにし、別の構造解析用コピーをアラインメントしてください。RMSF はアラインメントを行いません。ボックスが回転していない回転済みトラジェクトリに対して周期的コンタクトを実行しないでください。アラインメント、PCA、DSSP、水素結合、および母集団由来の自由エネルギーについては、解析リファレンスを参照してください。

def compute_rmsf(u, selection="protein and name CA", start_frame=0):
    """
    Compute per-residue RMSF (flexibility).

    Returns:
        residue_keys, rmsf_values; keys are (resindex, segid, resid, resname)
    """
    # Select atoms
    atoms = u.select_atoms(selection)

    if not len(atoms):
        raise ValueError("RMSF selection is empty")
    # A multi-atom selection returns mean atomic RMSF, not COM RMSF.
    R = rms.RMSF(atoms)
    R.run(start=start_frame)

    # Average by residue
    resids = []
    rmsf_per_res = []
    for res in u.select_atoms(selection).residues:
        res_atoms = res.atoms.intersection(atoms)
        if len(res_atoms) > 0:
            resids.append((res.ix, res.segid, res.resid, res.resname))
            # RMSF is indexed within the selected AtomGroup, not the Universe.
            rmsf_per_res.append(R.results.rmsf[atoms.resindices == res.ix].mean())

    return resids, np.array(rmsf_per_res)

4. タンパク質-リガンドコンタクト

def analyze_contacts(u, protein_sel="protein", ligand_sel="resname LIG",
                      radius=4.5, start_frame=0, periodic=True):
    """
    Track protein-ligand contacts over trajectory.

    Args:
        radius: Atom-pair cutoff in Angstroms; any pair defines a residue contact.
        periodic: Require box data and use minimum-image distances.
    Returns sets of (resindex, segid, resid, resname), one per analyzed frame.
    """
    from MDAnalysis.lib.distances import distance_array

    protein = u.select_atoms(protein_sel)
    ligand = u.select_atoms(ligand_sel)

    if not len(protein) or not len(ligand):
        raise ValueError("Protein and ligand selections must both be nonempty")
    if not np.isfinite(radius) or radius <= 0:
        raise ValueError("radius must be finite and positive")
    contact_frames = []
    for ts in u.trajectory[start_frame:]:
        if periodic and (ts.dimensions is None or not np.all(np.isfinite(ts.dimensions))
                         or np.any(ts.dimensions[:3] <= 0)):
            raise ValueError("Periodic contacts require valid frame box dimensions")
        # This includes hydrogen atoms unless the caller selects heavy atoms.
        distances = contacts.contact_matrix(
            distance_array(protein.positions, ligand.positions, box=ts.dimensions if periodic else None),
            radius,
        )
        contact_residues = set()
        for i in range(distances.shape[0]):
            if distances[i].any():
                res = protein[i].residue
                contact_residues.add((res.ix, res.segid, res.resid, res.resname))
        contact_frames.append(contact_residues)

    return contact_frames

力場と準備

科学的な系に対して、検証済みのタンパク質／リガンド／水／イオンの組み合わせを選択してください。 動作する例では AMBER14/TIP3P-FB を採用していますが、これは普遍的な推奨ではありません。 現在の OpenMM は amber19-all.xml（ff19SB、DNA OL21、RNA OL3、lipid21）も同梱しており、 さらに独自の charmm36_2024/water.xml を伴う charmm36_2024.xml も同梱しています。汎用の TIP3P と CHARMM 修飾 TIP3P は互換ではありません。IDP アンサンブルは特に タンパク質－水のバランスに敏感です。4サイト水は、上記の3サイトの例とは異なり、 追加の粒子と一致する Modeller の幾何学を必要とします。

AMBER のタンパク質 XML は任意のリガンドをパラメータ化しません。検証済みの GAFF2 または OpenFF の経路を、明示的な立体化学、プロトン化、結合次数、原子 マッピング、電荷法とともに使用してください。リガンドのみの Interchange はタンパク質－リガンド系ではありません。 保守的な PDBFixer と OpenFF の例、インストール要件、現在のリリース詳細については、 準備リファレンス を参照してください。

検証とプロベナンス

• まず最小化から始めてください。有限のエネルギーは有効な構造を確立しません。
• 段階的な NVT → NPT 平衡化 → 本番実行を使用してください。この例には拘束がありません。
• H を含む結合を拘束した 2 fs ステップは出発点の選択です。より大きな ステップ／HMR には、その系に対する積分器、安定性、観測量の検証が必要です。
• 周期的な距離にはフレームごとのボックスを使用してください。信頼できる結合トポロジーを用いて 構造観測量のために分子を完全なものにしてください。アンラッピングでは欠落した結合を推測できません。
• バーンイン、自己相関、有効サンプル、レプリカ、不確かさを評価してください。 短時間の安定した温度／密度のトレースはコンフォメーションサンプリングを確立しません。
• トポロジー／原子順序、System/Integrator XML、力場ファイル／バージョン、 シード、プラットフォーム／精度、パラメータ、解析の選択を保存してください。チェックポイントは プラットフォーム／ハードウェア／バージョンに固有です。XML 状態はより移植性が高いですが、 すべての内部乱数生成器の状態を保持するわけではありません。実際のセットアップの再起動をテストしてください。

追加リソース

• OpenMM ドキュメント: https://openmm.org/documentation.html
• MDAnalysis ユーザーガイド: https://docs.mdanalysis.org/
• GROMACS（代替 MD エンジン）: https://manual.gromacs.org/
• NAMD（代替）: https://www.ks.uiuc.edu/Research/namd/
• CHARMM-GUI（ウェブベースのシステムビルダー）: https://charmm-gui.org/
• AmberTools（無料の Amber ツール）: https://ambermd.org/AmberTools.php
• OpenMM 論文: Eastman P et al. (2017) PLOS Computational Biology. PMID: 28278240
• MDAnalysis 論文: Michaud-Agrawal N et al. (2011) J Computational Chemistry. PMID: 21500218