分子動力學

分子動力學 是做什麼的？

概述

分子動力學模擬通過牛頓運動方程，以計算方式模擬分子系統隨時間的演化。此技能涵蓋兩個互補的工具：

• OpenMM (https://openmm.org/)：高效能 MD 模擬引擎，支援 GPU、Python API 和彈性的力場設定
• 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, topology
    """
    # 載入 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"系統：{modeller.topology.getNumAtoms()} 個原子，"
          f"{modeller.topology.getNumResidues()} 個殘基")

    # 創建系統
    system = forcefield.createSystem(
        modeller.topology,
        nonbondedMethod=PME,         # 粒子網格 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 系統
        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"初始能量：{state.getPotentialEnergy()}")

    # 最小化
    simulation.minimizeEnergy(
        tolerance=tolerance*kilojoules_per_mole/nanometer,
        maxIterations=max_iterations
    )

    state = simulation.context.getState(getEnergy=True, getPositions=True)
    print(f"最小化後能量：{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 模擬物件（最小化後）
        n_steps: MD 步數 (50000 × 2 fs = 100 ps)
        temperature: 溫度 (Kelvin)
        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"執行 NVT 平衡：{n_steps} 步 ({n_steps*2/1000:.1f} ps)")
    simulation.step(n_steps)
    print("NVT 平衡完成")

    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 × 2 fs = 1 ns)
        temperature: 溫度 (Kelvin)
        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"執行 NPT 生產：{n_steps} 步 ({n_steps*2/1000000:.2f} ns)")
    simulation.step(n_steps)
    print("生產 MD 完成")
    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} 個原子, {u.trajectory.n_frames} 幀")
    print(f"時間範圍: 0 到 {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:
        (時間, 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: 接觸距離截止值 (Angstroms)
    """
    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 問題：缺失殘基、原子、添加 H、標準化。"""
    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-GUI（系統準備工具）：https://charmm-gui.org/