BioPlayground

🧬
목록으로

분자동역학 입문: GROMACS·OpenMM으로 단백질이 움직이는 걸 본다

정적인 구조 예측을 넘어, 원자가 시간에 따라 어떻게 움직이는지를 뉴턴 방정식으로 직접 풀어내는 분자동역학(MD)의 원리와 OpenMM 실습을 다룹니다.

심화
|
25분
|
검증 완료 (2026-07-29)
molecular dynamicsGROMACSOpenMMforce field
진행률0/120 (0%)

F09 다음 질문 — 도킹은 정지 사진 한 장이다

F09에서 AutoDock Vina로 리간드가 결합 부위에 어떤 자세로 붙는지 계산했습니다. 그런데 그 결과는 정확히 말하면 정지 사진 한 장입니다. 실제 단백질과 리간드는 체온에서 끊임없이 흔들리고, 결합 후에도 완전히 고정되지 않으며, 어떤 결합은 순간적으로 형성됐다 풀리기를 반복합니다. F02～F08에서 다룬 구조 예측·도킹 모델들이 공통으로 답하지 못하는 질문이 여기 있습니다 — 이 구조가 시간에 따라 어떻게 움직이는가.

**분자동역학(Molecular Dynamics, MD)**은 이 질문에 물리학의 가장 기본적인 도구, 뉴턴의 운동방정식으로 직접 답합니다. GROMACS와 OpenMM은 이 계산을 실무에서 돌리는 대표적인 두 엔진입니다.

원리 — 뉴턴 방정식을 원자 하나하나에 적용한다

힘장(Force Field): 원자 사이의 힘을 정의하는 규칙

MD의 출발점은 모든 원자 쌍 사이에 작용하는 힘을 계산하는 것입니다. 이 힘을 정의하는 규칙의 집합을 **힘장(force field)**이라 부릅니다. F09에서 다룬 도킹의 반데르발스·정전기 항과 비슷한 구성이지만, MD 힘장은 여기에 결합 길이·각도·이면각(dihedral)까지 포함한 결합 내부 항까지 함께 다룹니다.

Utotal=∑bondskb(r−r0)2⏟결합 신축+∑angleskθ(θ−θ0)2⏟결합각 굽힘+∑dihedralskϕ[1+cos⁡(nϕ−δ)]⏟이면각 회전+∑i<j[Aijrij12−Bijrij6+qiqj4πϵ0rij]⏟비결합 상호작용U_{\text{total}} = \underbrace{\sum_{\text{bonds}} k_b(r-r_0)^2}_{\text{결합 신축}} + \underbrace{\sum_{\text{angles}} k_\theta(\theta-\theta_0)^2}_{\text{결합각 굽힘}} + \underbrace{\sum_{\text{dihedrals}} k_\phi[1+\cos(n\phi-\delta)]}_{\text{이면각 회전}} + \underbrace{\sum_{i<j}\left[\frac{A_{ij}}{r_{ij}^{12}} - \frac{B_{ij}}{r_{ij}^{6}} + \frac{q_iq_j}{4\pi\epsilon_0 r_{ij}}\right]}_{\text{비결합 상호작용}}

앞의 세 항은 원자가 화학결합으로 이어진 이웃끼리의 진동·회전을 스프링 모델로 근사한 것이고, 마지막 항은 F09의 Lennard-Jones + 쿨롱 항과 같은 비결합 상호작용입니다. AMBER, CHARMM 같은 이름을 들어봤다면, 그게 바로 이 항들의 파라미터 kb,kθ,Aij,Bijk_b, k_\theta, A_{ij}, B_{ij} 등을 실험·양자화학 계산으로 정밀하게 피팅해둔 구체적인 힘장 세트입니다.

뉴턴 방정식과 수치 적분

이 퍼텐셜 에너지 UU로부터 각 원자에 작용하는 힘은 에너지의 음의 기울기입니다.

Fi=−∇iU,mid2ridt2=FiF_i = -\nabla_i U, \qquad m_i \frac{d^2 r_i}{dt^2} = F_i

이건 순수한 뉴턴의 제2법칙입니다. 문제는 원자 수천～수만 개가 서로 얽힌 이 미분방정식 시스템을 해석적으로 풀 수 없다는 점입니다. 그래서 시간을 아주 작은 간격(Δt\Delta t, 보통 1～2펨토초)으로 쪼개, 각 스텝마다 위치와 속도를 갱신하는 수치 적분을 씁니다. 가장 널리 쓰이는 것이 leapfrog 또는 velocity Verlet 적분법입니다.

v(t+Δt2)=v(t−Δt2)+F(t)mΔtv\left(t+\frac{\Delta t}{2}\right) = v\left(t-\frac{\Delta t}{2}\right) + \frac{F(t)}{m}\Delta t

r(t+Δt)=r(t)+v(t+Δt2)Δtr(t+\Delta t) = r(t) + v\left(t+\frac{\Delta t}{2}\right)\Delta t

속도를 반 스텝 어긋나게(leapfrog, "개구리 뛰기") 갱신하는 이 방식은 에너지 보존 성질이 좋아 장시간 시뮬레이션에서 오차가 누적되지 않는 장점이 있습니다.

손 계산 예제: 왜 타임스텝이 펨토초 단위여야 하는가

수소 원자가 관여하는 결합 진동의 주기는 대략 10 펨토초(10−1410^{-14}초) 안팎입니다. 수치 적분이 안정적으로 작동하려면, 타임스텝 Δt\Delta t가 이 가장 빠른 진동 주기보다 충분히 작아야 합니다. 대략적인 안정성 기준으로 진동 주기의 약 1/10 이하를 타임스텝으로 잡는다고 하면,

Δt≲10 fs10=1 fs\Delta t \lesssim \frac{10\,\text{fs}}{10} = 1\,\text{fs}

즉 타임스텝을 1～2펨토초로 잡는 이유는 임의의 관행이 아니라, 시스템에서 가장 빠르게 진동하는 결합(주로 수소가 관여하는 결합)을 놓치지 않기 위한 물리적 제약입니다. 만약 1나노초(10−910^{-9}초) 분량의 궤적을 얻고 싶다면 필요한 스텝 수는

1 ns1 fs=10−910−15=106 steps\frac{1\,\text{ns}}{1\,\text{fs}} = \frac{10^{-9}}{10^{-15}} = 10^{6}\,\text{steps}

100만 스텝입니다. 이 편의 제목에 "100ns 궤적"이 언급된 것처럼, 실무에서 의미 있는 시간 스케일(수십～수백 나노초)에 도달하려면 스텝 수가 억 단위로 늘어나므로, MD가 왜 그렇게 계산 집약적인지 이 손 계산에서 바로 체감할 수 있습니다.

실습: OpenMM으로 짧은 펩타이드 시뮬레이션

python
# Colab T4에서 실행
!pip install -q openmm


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


# 예시: Colab 세션에 실제로 내려받은 작은 결정구조 crambin(1CRN)
from urllib.request import urlretrieve
urlretrieve("https://files.rcsb.org/download/1CRN.pdb", "1CRN.pdb")
pdb = PDBFile("1CRN.pdb")
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')


# 결정구조에 빠진 수소를 force-field template에 맞춰 추가합니다.
modeller = Modeller(pdb.topology, pdb.positions)
modeller.addHydrogens(forcefield, pH=7.0)


# 교육용 진공 시뮬레이션입니다. 실제 분석에는 용매·이온·PME 설정이 필요합니다.
system = forcefield.createSystem(modeller.topology, nonbondedMethod=NoCutoff, constraints=HBonds)


integrator = LangevinMiddleIntegrator(
    300 * kelvin,        # 목표 온도
    1 / picosecond,       # 마찰 계수 (열욕과의 결합 강도)
    2 * femtoseconds       # 타임스텝 — 앞의 손 계산과 정확히 같은 규모
)


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


# 1) 에너지 최소화 — 초기 구조의 비정상적 원자 겹침(steric clash)을 먼저 제거
simulation.minimizeEnergy()


# 2) 짧은 production run — 1만 스텝 = 20 피코초 (교육용 축소 규모)
simulation.reporters.append(
    StateDataReporter(stdout, 1000, step=True, potentialEnergy=True, temperature=True)
)
simulation.step(10000)

constraints=HBonds는 수소가 관여하는 결합의 진동을 아예 고정(constraint)해버리는 트릭입니다. 앞의 손 계산에서 본 것처럼 가장 빠른 진동이 타임스텝을 제약하는데, 이 진동 자체를 걸어 잠그면 타임스텝을 2펨토초까지 늘릴 수 있어 실무에서 널리 씁니다.

CS 매핑

  • 수치 적분(Numerical Integration): 연속 미분방정식을 이산 시간 스텝으로 근사해 푸는 leapfrog/Verlet 적분법은, 물리 엔진·그래픽스의 강체 시뮬레이션에서 쓰는 수치 적분기와 원리적으로 동일합니다.
  • 상태 공간 탐색(State Space Exploration): 시뮬레이션이 매 스텝 시스템의 위치·속도(상태)를 갱신하며 궤적을 따라가는 과정은, 강화학습에서 에이전트가 상태 공간을 순차적으로 탐색하는 것과 구조적으로 유사합니다.
  • 제약 조건을 통한 연산 절약(Constraint-based Speedup): 빠른 진동을 constraint로 고정해 타임스텝을 늘리는 트릭은, 계산 그래프에서 세밀한 부분을 근사·고정해 전체 연산을 절약하는 일반적인 최적화 패턴입니다.

자주 만나는 결함

  • 에너지 최소화 없이 바로 시뮬레이션 시작: 초기 구조에 원자 겹침이 남아 있으면 첫 스텝부터 힘이 비정상적으로 커져 시스템이 발산(explosion)합니다. minimizeEnergy()는 생략하면 안 되는 필수 단계입니다.
  • 평형화(equilibration) 없이 production run 결과를 바로 분석: 초기 구조는 목표 온도·압력에 아직 적응하지 못한 상태입니다. 충분한 평형화 구간을 거치지 않은 궤적 초반부를 분석에 포함하면 인위적인 편향이 섞입니다.
  • 짧은 시뮬레이션에서 통계적으로 유의미한 결론 도출: 이 실습의 20피코초는 순전히 교육용 축소 규모입니다. 실제 결합 안정성이나 구조 변화를 논하려면 최소 수십～수백 나노초, 경우에 따라 여러 반복(replica)이 필요합니다.
  • 진공 예제를 실제 생체 환경으로 오해: 위 코드는 API와 적분 흐름을 보여주는 최소 예제입니다. 논문 수준 계산에는 명시적 물·이온, 주기 경계조건, 장거리 정전기(PME), NVT/NPT 평형화, force-field/ligand parameter 검증과 독립 반복이 필요합니다.

더 깊게 파고 싶다면

본문은 BPD 연구진이 직접 재구성한 서술입니다. 원 논문과 공식 자료로 심화해봅시다.

  • GROMACS 공식 문서·튜토리얼: manual.gromacs.org — 힘장 선택, 시스템 준비 전체 워크플로.
  • OpenMM 공식 문서: openmm.org — Python API와 GPU 가속 설정.
  • AMBER 힘장 논문: Maier et al. (2015), ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters, J. Chem. Theory Comput.

F1.2(오픈 구조예측·도킹·MD) 5편을 마쳤습니다. 다음 5편(F11～F15)에서는 다시 구조가 아니라 서열 자체로 돌아가, 단백질 언어모델이 임베딩·생성·설계에 어떻게 쓰이는지 다룹니다.