molecular-dynamics

分子動力学

概要

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

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

インストール方法：

conda install -c conda-forge openmm mdanalysis nglview
# または
uv pip install openmm mdanalysis

このスキルを使用すべき場合

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

• タンパク質安定性解析：変異はタンパク質の動的挙動にどのように影響するか？
• 薬剤結合シミュレーション：リガンドの結合モードと滞留時間を特徴づける
• 構造サンプリング：タンパク質の柔軟性と構造変化を探る
• タンパク質間相互作用：界面ダイナミクスと結合エネルギーをモデル化する
• RMSD/RMSF解析：基準構造からの構造的ゆらぎを定量化する
• 自由エネルギー推定：結合自由エネルギーまたはコンフォメーション自由エネルギーを計算する
• 膜シミュレーション：脂質二重層中のタンパク質をモデル化する
• 天然変性タンパク質：IDRのコンフォメーションアンサンブルを研究する

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

1. システム準備

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

def prepare_system_from_pdb(pdb_file, forcefield_name="amber14-all.xml",
                              water_model="amber14/tip3pfb.xml"):
    """
    PDBファイルからOpenMMシステムを準備する。

    Args:
        pdb_file: クリーニング済みPDBファイルへのパス（生のPDBファイルにはPDBFixerを使用）
        forcefield_name: 力場XMLファイル
        water_model: 水モデルXMLファイル

    Returns:
        pdb, forcefield, system, topography
    """
    # PDBの読み込み
    pdb = PDBFile(pdb_file)

    # 力場の読み込み
    forcefield = ForceField(forcefield_name, water_model)

    # 水素原子の追加と溶媒和
    modeller = Modeller(pdb.topology, pdb.positions)
    modeller.addHydrogens(forcefield)

    # 溶媒ボックスの追加（10 Åのパディング、150 mM NaCl）
    modeller.addSolvent(
        forcefield,
        model='tip3p',
        padding=10*angstroms,
        ionicStrength=0.15*molar
    )

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

    # システムの作成
    system = forcefield.createSystem(
        modeller.topology,
        nonbondedMethod=PME,         # 長距離静電相互作用にはParticle Mesh Ewald
        nonbondedCutoff=1.0*nanometer,
        constraints=HBonds,           # 水素結合を拘束（2 fsのタイムステップを可能にする）
        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):
    """
    系から立体障害を取り除くためにエネルギー最小化を行う。

    Args:
        modeller: トポロジーと位置を持つModellerオブジェクト
        system: OpenMM System
        output_pdb: 最小化された構造の保存先パス
        max_iterations: 最小化ステップの最大数
        tolerance: 収束基準（kJ/mol/nm）

    Returns:
        最小化された位置を持つシミュレーションオブジェクト
    """
    # 積分器の設定（最小化には関係しない）
    integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)

    # シミュレーションの作成
    # 可能であればGPUを使用（CUDAまたはOpenCL）、CPUにフォールバック
    try:
        platform = Platform.getPlatformByName('CUDA')
        properties = {'DeviceIndex': '0', 'Precision': 'mixed'}
    except Exception:
        try:
            platform = Platform.getPlatformByName('OpenCL')
            properties = {}
        except Exception:
            platform = Platform.getPlatformByName('CPU')
            properties = {}

    simulation = Simulation(
        modeller.topology, system, integrator,
        platform, properties
    )
    simulation.context.setPositions(modeller.positions)

    # 初期エネルギーの確認
    state = simulation.context.getState(getEnergy=True)
    print(f"Initial energy: {state.getPotentialEnergy()}")

    # 最小化
    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()}")

    # 最小化された構造を保存
    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平衡化：一定のN, V, T。
    速度を目標温度に平衡化する。

    Args:
        simulation: OpenMM Simulation（最小化後）
        n_steps: MDステップ数（50000 × 2fs = 100 ps）
        temperature: 温度（ケルビン）
        report_interval: データレポート間のステップ数
        output_prefix: トラジェクトリとログのファイル接頭辞
    """
    # NVT中は主鎖に位置拘束を追加
    # （オプション：重原子を拘束）

    # 温度設定
    simulation.context.setVelocitiesToTemperature(temperature*kelvin)

    # レポーター追加
    simulation.reporters = []

    # ログファイル
    simulation.reporters.append(
        StateDataReporter(
            f"{output_prefix}_log.txt",
            report_interval,
            step=True,
            potentialEnergy=True,
            kineticEnergy=True,
            temperature=True,
            volume=True,
            speed=True
        )
    )

    # DCDトラジェクトリ（コンパクトなバイナリ形式）
    simulation.reporters.append(
        DCDReporter(f"{output_prefix}_traj.dcd", report_interval)
    )

    print(f"Running NVT equilibration: {n_steps} steps ({n_steps*2/1000:.1f} ps)")
    simulation.step(n_steps)
    print("NVT equilibration complete")

    return simulation

4. NPT平衡化とプロダクションラン

def run_npt_production(simulation, n_steps=500000, temperature=300, pressure=1.0,
                        report_interval=5000, output_prefix="npt"):
    """
    NPTプロダクションラン：一定のN, P, T。

    Args:
        n_steps: プロダクションステップ数（500000 × 2fs = 1 ns）
        temperature: 温度（ケルビン）
        pressure: 圧力（bar）
        report_interval: レポート間のステップ数
    """
    # 圧力制御のためのモンテカルロバロスタットを追加
    system = simulation.context.getSystem()
    system.addForce(MonteCarloBarostat(pressure*bar, temperature*kelvin, 25))
    simulation.context.reinitialize(preserveState=True)

    # レポーターを更新
    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)
    )

    # チェックポイントを保存
    simulation.reporters.append(
        CheckpointReporter(f"{output_prefix}_checkpoint.chk", 50000)
    )

    print(f"Running NPT production: {n_steps} steps ({n_steps*2/1000000:.2f} ns)")
    simulation.step(n_steps)
    print("Production MD complete")
    return simulation

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

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

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):
    """
    MDAnalysisでMDトラジェクトリを読み込む。

    Args:
        topology_file: PDB、PSF、またはその他のトポロジーファイル
        trajectory_file: DCD、XTC、TRR、またはその他のトラジェクトリ
    """
    u = mda.Universe(topology_file, trajectory_file)
    print(f"Universe: {u.atoms.n_atoms} atoms, {u.trajectory.n_frames} frames")
    print(f"Time range: 0 to {u.trajectory.totaltime:.0f} ps")
    return u

2. RMSD解析

def compute_rmsd(u, selection="backbone", reference_frame=0):
    """
    選択された原子のRMSDを参照フレームに対して計算する。

    Args:
        u: MDAnalysis Universe
        selection: 原子選択文字列（MDAnalysis構文）
        reference_frame: 基準構造のフレームインデックス

    Returns:
        (time, rmsd) 値のnumpy配列
    """
    # トラジェクトリを整列してRMSDを最小化
    aligner = align.AlignTraj(u, u, select=selection, in_memory=True)
    aligner.run()

    # RMSDの計算
    R = rms.RMSD(u, select=selection, ref_frame=reference_frame)
    R.run()

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

def plot_rmsd(rmsd_data, title="RMSD over time", output_file="rmsd.png"):
    """シミュレーション時間に対するRMSDをプロットする。"""
    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解析（残基ごとの柔軟性）

def compute_rmsf(u, selection="backbone", start_frame=0):
    """
    残基ごとのRMSF（柔軟性）を計算する。

    Returns:
        resids, rmsf_values 配列
    """
    # 原子の選択
    atoms = u.select_atoms(selection)

    # RMSFの計算
    R = rms.RMSF(atoms)
    R.run(start=start_frame)

    # 残基ごとに平均化
    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.resid)
            rmsf_per_res.append(R.results.rmsf[res_atoms.indices].mean())

    return np.array(resids), np.array(rmsf_per_res)

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

def analyze_contacts(u, protein_sel="protein", ligand_sel="resname LIG",
                      radius=4.5, start_frame=0):
    """
    トラジェクトリ全体でのタンパク質-リガンドコンタクトを追跡する。

    Args:
        radius: コンタクト距離カットオフ（Å）
    """
    protein = u.select_atoms(protein_sel)
    ligand = u.select_atoms(ligand_sel)

    contact_frames = []
    for ts in u.trajectory[start_frame:]:
        # リガンドからの半径内にあるタンパク質原子を見つける
        distances = contacts.contact_matrix(
            protein.positions, ligand.positions, radius
        )
        contact_residues = set()
        for i in range(distances.shape[0]):
            if distances[i].any():
                contact_residues.add(protein.atoms[i].resid)
        contact_frames.append(contact_residues)

    return contact_frames

力場選択ガイド

システム	推奨力場	水モデル
標準的なタンパク質	AMBER14（amber14-all.xml）	TIP3P-FB
タンパク質 + 小分子	AMBER14 + GAFF2	TIP3P-FB
膜タンパク質	CHARMM36m	TIP3P
核酸	AMBER99-bsc1 または AMBER14	TIP3P
天然変性タンパク質	ff19SB または CHARMM36m	TIP3P

システム準備ツール

PDBFixer（生のPDBファイル用）

from pdbfixer import PDBFixer
from openmm.app import PDBFile

def fix_pdb(input_pdb, output_pdb, ph=7.0):
    """一般的なPDB問題を修正：欠落残基、原子、水素追加、標準化。"""
    fixer = PDBFixer(filename=input_pdb)
    fixer.findMissingResidues()
    fixer.findNonstandardResidues()
    fixer.replaceNonstandardResidues()
    fixer.removeHeterogens(True)    # 水/リガンドを除去
    fixer.findMissingAtoms()
    fixer.addMissingAtoms()
    fixer.addMissingHydrogens(ph)

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

    return output_pdb

小分子用GAFF2（OpenFF Toolkit経由）

# リガンドパラメータ化には、OpenFF ToolkitまたはACPYPEを使用
# 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"):
    """小分子用のGAFF2/OpenFFパラメータを生成する。"""
    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

ベストプラクティス

• MDの前に必ず最小化を行う：生のPDB構造には立体障害がある
• プロダクションの前に平衡化する：NVT（50–100 ps）→ NPT（100–500 ps）→ プロダクション
• GPUを使用する：GPU（CUDA/OpenCL）でのシミュレーションは10～100倍高速
• HBonds拘束による2 fsのタイムステップ：標準。HMR（水素質量再分配）では4 fsを使用
• 平衡化されたトラジェクトリのみを解析する：最初の20～50%は平衡化として破棄する
• チェックポイントを保存する：MD実行は失敗することがある。チェックポイントにより再開可能
• 周期境界条件：溶媒和系に必要
• 静電相互作用にはPME：荷電系ではカットオフ法よりも正確

追加リソース

• 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-G…