분자동역학

분자동역학은(는) 무엇을 하나요?

개요

분자동역학(MD) 시뮬레이션은 뉴턴의 운동 방정식을 적분하여 분자 시스템의 시간 진화를 계산적으로 모델링합니다. 이 스킬은 두 가지 보완적인 도구를 다룹니다:

• 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, 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,         # 장거리 정전기학을 위한 입자 메시 에발트
        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"초기 에너지: {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 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"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 × 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"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  # 열: 프레임, 시간, RMSD
    return rmsd_data

def plot_rmsd(rmsd_data, title="시간에 따른 RMSD", 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("시간 (ns)")
    ax.set_ylabel("RMSD (Å)")
    ax.set_title(title)
    ax.axhline(rmsd_data[:, 2].mean(), color='r', linestyle='--',
               label=f'평균: {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 문제 수정: 누락된 잔기, 원자, 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 툴킷 사용)

# 리간드 파라미터화를 위해 OpenFF 툴킷 또는 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에서 10–100배 더 빠릅니다 (CUDA/OpenCL).
• 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://www.charmm-gui.org/