분자동역학

1. 개요

분자동역학 Molecular Dynamics
약칭	MD
분야	계산화학 × 통계역학 × 생물물리
지배 법칙	뉴턴 운동방정식 F = ma
시간 적분	속도 베를레(심플렉틱)
시간 스텝	1 fs (구속 쓰면 2 fs)
장거리 쿨롱	에발트 합산 / PME, O(N log N)
대표 소프트웨어	LAMMPS, GROMACS

결국 원자 하나하나에 F=ma 때려박고 시간을 흘려보내는 것, 그 이상도 이하도 아니다.

분자동역학(Molecular Dynamics, MD)은 원자와 분자로 이루어진 계에서 각 입자에 뉴턴의 운동방정식을 적용하고 이를 시간에 대해 수치적으로 적분함으로써, 계가 시간에 따라 어떻게 변하는지를 추적하는 시뮬레이션 기법이다. 양자역학적 슈뢰딩거 방정식을 푸는 대신, 원자를 고전적인 점입자로 취급하고 그 궤적을 좇는다는 점에서 계산량이 극적으로 줄어든다.¹

단백질의 접힘(folding), 신약 후보와 단백질의 결합, 고분자의 물성, 유리의 구조, 나노 재료의 열전달 — 원자 스케일에서 “시간에 따라 무슨 일이 벌어지는가”를 알고 싶을 때 MD가 등장한다. 전산유체역학이 연속체를 격자로 쪼갠다면, MD는 물질을 아예 원자 알갱이 단위로 쪼개서 하나하나 굴린다. 현대적 형태의 출발점은 1967년 베를레(Verlet)가 레너드-존스 퍼텐셜 입자 864개를 적분해 아르곤의 열역학 물성을 재현한 논문으로, 오늘날 쓰는 적분기와 이웃 목록(Verlet list)이 둘 다 거기서 나왔다.

2. 지배 방정식과 힘장

MD의 심장은 지극히 고전적이다. $i$번째 원자의 위치 $\mathbf{r}_i$는 뉴턴 제2법칙을 따른다:

$$m_i \frac{d^2 \mathbf{r}_i}{dt^2} = \mathbf{F}_i = -\nabla_i U(\mathbf{r}_1, \dots, \mathbf{r}_N)$$

핵심은 퍼텐셜 에너지 함수 $U$, 이른바 힘장(force field)이다. 이게 원자들 사이의 상호작용을 규정하며, MD 결과의 신뢰도는 전적으로 힘장의 품질에 달려 있다. 대표적으로 비결합 상호작용을 기술하는 레너드-존스 퍼텐셜(Lennard-Jones potential)은:

$$U_{\text{LJ}}(r) = 4\varepsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right]$$

가까우면 강하게 밀치고($r^{-12}$ 반발), 적당히 멀면 살짝 끌어당기는($r^{-6}$ 인력) 형태다. 최소점은 $r = 2^{1/6}\sigma = 1.1225\,\sigma$ 에서 $U = -\varepsilon$ 이다. 여기에 결합 신축, 결합각, 이면각, 정전기 상호작용 항들이 더해져 AMBER, CHARMM, OPLS 같은 완성된 힘장을 이룬다.²

비결합 항은 쌍이 $N(N-1)/2$ 개라 그냥 다 더하면 $O(N^2)$ 다. LJ 항은 짧은 거리에서 급히 죽으므로 절단 반경 $r_c$(보통 $2.5\sigma$ 또는 10~12 Å)로 잘라 이웃 목록과 셀 분할로 $O(N)$ 에 처리한다. 반면 쿨롱 항은 $1/r$ 로 느리게 죽어 절단이 허용되지 않고, 에발트 합산을 FFT로 가속한 PME(particle mesh Ewald)로 $O(N\log N)$ 에 푼다. 이 두 가지가 실제 MD 런타임의 대부분을 먹는다.

3. 시간 적분

운동방정식을 컴퓨터로 푸는 것은 결국 이산적인 시간 스텝마다 위치와 속도를 업데이트하는 일이다. MD에서 압도적으로 널리 쓰이는 것은 베를레 적분(Verlet integration) 계열, 특히 속도 베를레(velocity Verlet) 알고리즘이다:

$$\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\,\Delta t + \frac{1}{2}\mathbf{a}(t)\,\Delta t^2$$ $$\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{1}{2}\left[\mathbf{a}(t) + \mathbf{a}(t + \Delta t)\right]\Delta t$$

베를레류가 사랑받는 이유는 시간 가역성과 심플렉틱(symplectic) 성질 때문이다. 국소 정확도는 겨우 2차인데도 전체 에너지가 장기적으로 표류하지 않고 진동만 하므로, 정교한 룽게-쿠타법보다 오히려 장시간 안정성이 좋다. 화려함보다 보존성이 이기는 드문 사례이고, 아래 § 검산에서 이걸 숫자로 확인한다.

4. 시간 스텝의 저주

MD의 근본적 한계는 시간에 있다. 가장 빠른 운동이 수소 원자의 결합 신축(O–H, C–H)인데, 진동수가 약 30003700 cm⁻¹ 라 주기가 **911 fs** 밖에 안 된다. 이걸 제대로 잡으려면 시간 스텝 $\Delta t$ 를 주기의 1/10 수준인 1 fs($10^{-15}$초)로 잘게 써야 한다.³

그런데 우리가 보고 싶은 현상들의 시간 스케일은:

현상	시간 스케일	필요 스텝 수 (1 fs)
결합 진동	~10 fs	~10
물 분자 재배열	~1 ps	~1,000
작은 단백질 접힘	~1 μs	~10억
큰 단백질/막 과정	~1 ms	~1조

1 밀리초짜리 현상을 보려면 스텝을 1조 번 밟아야 한다는 뜻이다. 이 “시간 스케일 격차”가 MD의 아킬레스건이며, 그래서 전용 슈퍼컴퓨터(D. E. Shaw의 Anton 등)나 향상된 샘플링 기법(자유에너지 섭동, 복제 교환 등)이 동원된다. 공간이 아니라 시간을 병렬 컴퓨팅과 GPU 컴퓨팅에 갈아 넣는 셈이다 — 다만 시간 방향은 본질적으로 순차라서, 병렬화로 줄어드는 것은 “스텝 하나의 벽시계 시간”뿐이고 스텝 수 자체는 줄지 않는다. 이것이 MD 병렬화의 천장이다.

5. 검산

시간 스텝의 근거. 진동수 $\tilde\nu$ [cm⁻¹]의 주기는 $T = 1/(c\tilde\nu)$ 다. 직접 계산하면:

진동 모드	$\tilde\nu$ (cm⁻¹)	주기
O–H 신축	3700	9.0 fs
C–H 신축	3000	11.1 fs
결합각 굽힘	1700	19.6 fs
중원자 신축	1000	33.4 fs

$\Delta t = 1$ fs 면 가장 빠른 진동을 주기당 9점으로 샘플링한다 — 심플렉틱 적분기가 안정하게 버티는 하한 근처다. SHAKE/LINCS로 수소 결합 길이를 구속해 이 모드를 제거하면 다음으로 빠른 것이 ~20 fs 주기이므로 $\Delta t = 2$ fs 가 정당화된다. 2 fs 로 바꿔도 1 ms 를 보려면 여전히 $5\times10^{11}$ 스텝이고, 1조에서 5천억으로 준 것이 전부다. 구속은 정확도 조정이 아니라 2배 할인이라는 뜻이다.

심플렉틱이 왜 이기는가. 조화진동자 $\ddot x = -x$ ($E_0 = 0.5$)를 $\Delta t = 0.1$ 로 적분했을 때 에너지 상대오차:

스텝 수	속도 베를레	RK4	양함수 오일러
$10^3$	$-5.5\times10^{-4}$	$-1.4\times10^{-5}$	$+2.1\times10^{4}$
$10^5$	$-2.4\times10^{-3}$	$-1.4\times10^{-3}$	발산
$10^7$	$-1.9\times10^{-3}$	$-13\%$	NaN

짧게 보면 4차 정확도인 RK4가 베를레보다 40배 정확하다. 그런데 $10^7$ 스텝(약 16만 주기)에서는 RK4의 에너지가 13% 흘러내린 반면, 베를레는 $10^{-3}$ 수준에서 더 커지지 않고 진동만 한다(오히려 $10^5$ 스텝보다 작다). MD가 수억 스텝을 밟는 도구라는 점을 떠올리면, 어느 쪽을 골라야 하는지는 토론의 여지가 없다. 양함수 오일러는 $10^3$ 스텝 만에 에너지를 4자리 불려 터뜨린다.

온도는 속도에서 나온다. 등분배 정리로 $\langle E_k\rangle = \frac{3}{2}Nk_BT$ 이므로 온도는 측정량이 아니라 속도의 통계다. 아르곤($m = 39.95$ u) 300 K 에서 $v_{\text{rms}} = \sqrt{3k_BT/m} = 433$ m/s — 1 fs 동안 0.43 pm 를 간다. 셀 크기(Å 단위)에 비해 턱없이 작은 이 비율이 MD가 스텝을 수십억 번 밟아야 하는 물리적 이유다. 거꾸로 LJ 자연단위 시간 $\tau = \sigma\sqrt{m/\varepsilon}$ 을 아르곤 값($\sigma = 3.405$ Å, $\varepsilon/k_B = 119.8$ K)으로 환산하면 $\tau = 2.16$ ps 이고, 보통 $\Delta t = 0.005\tau \approx 10.8$ fs 를 쓴다. 수소가 없는 계에서는 스텝을 열 배 키울 수 있다는 뜻이다.

절단 반경의 대가. $r_c = 2.5\sigma$ 에서 LJ 퍼텐셜 값은 $-1.63\times10^{-2}\varepsilon$ 로 최소점의 1.6%다. $3.0\sigma$ 에서 $-5.5\times10^{-3}$, $3.5\sigma$ 에서 $-2.2\times10^{-3}$. 그냥 자르면 이 값만큼의 계단이 힘에 생겨 에너지가 표류하므로, 실무에서는 이동(shift)·스위칭이나 장거리 꼬리 보정을 반드시 켠다.

6. 앙상블과 주기경계조건

현실적인 시뮬레이션을 위해 몇 가지 장치가 필수적으로 붙는다.

• 앙상블 — 어떤 열역학 변수를 고정할지에 따라 나뉜다. NVE(입자수·부피·에너지 고정)는 순수 뉴턴 역학 그대로다. 실제 실험 조건에 맞추려면 온도를 고정한 NVT, 온도와 압력을 고정한 NPT를 쓴다.
• 온도조절기(thermostat) — NVT/NPT에서 온도를 목표값에 붙들어두는 알고리즘. 노제-후버(Nosé–Hoover), 랑주뱅(Langevin), 베렌센(Berendsen) 등이 있다. 온도란 결국 원자들의 평균 운동에너지이므로, 조절기는 속도를 슬쩍슬쩍 스케일링해 온도를 맞춘다. 단 베렌센류는 정준분포를 정확히 재현하지 못해(이른바 “날아가는 얼음덩이” 문제) 평형 샘플링에는 노제-후버나 확률적 조절기를 쓰는 것이 맞다.
• 주기경계조건(periodic boundary conditions, PBC) — 유한한 계산 상자로 무한한 벌크를 흉내 내는 트릭. 상자를 벗어난 원자는 반대편으로 다시 들어온다. 표면 효과 없이 대량 물질을 모사할 수 있게 해주는 필수 요소지만, 상자가 작으면 분자가 자기 복제본과 상호작용해 인공적인 유한크기 효과가 생긴다.

7. 몬테카를로와의 대비

MD와 자주 짝지어 언급되는 것이 몬테카를로 방법(Monte Carlo, MC)이다. 둘 다 통계역학적 평형 성질을 뽑아내는 데 쓰이지만 철학이 다르다.

MD는 시간을 따라 “굴린다”. MC는 시간을 버리고 무작위로 “찍는다”.

MD는 실제 물리적 궤적을 따라가므로 동역학적 정보(확산계수, 점성, 시간 상관함수 등)를 줄 수 있다. 반면 MC는 확률 분포를 무작위 표본추출로 훑기 때문에 시간 개념이 없고 오직 평형 성질만 준다. 대신 MC는 에너지 장벽을 껑충 뛰어넘는 이동이 가능해, MD가 갇혀버리는 상태 공간을 더 효율적으로 탐색하기도 한다. 실전에서는 둘을 섞은 하이브리드 기법(하이브리드 MC)도 흔하다.

8. 주요 소프트웨어

• LAMMPS — 재료·물성 분야의 오픈소스 강자. 공간 분할(spatial decomposition) 병렬화로 수십억 원자 규모까지 확장한다(Plimpton 1995가 그 설계 문서다).
• GROMACS — 생체분자(단백질·막·물) 시뮬레이션의 표준. SIMD 커널과 GPU 오프로드 최적화가 극악하기로 유명하다.
• NAMD, AMBER, OpenMM — 각각 대규모 병렬, 힘장 생태계, GPU/파이썬 친화성으로 특화된 코드들.

9. 관련 문서

• 베를레 적분 · 심플렉틱 적분기
• 레너드-존스 퍼텐셜 · 에발트 합산
• 랑주뱅 동역학 · 주기경계조건
• 몬테카를로 방법 · 자유에너지 섭동
• 슈뢰딩거 방정식 · 밀도범함수이론
• 룽게-쿠타법 · SPH
• LAMMPS · GROMACS
• 병렬 컴퓨팅 · GPU 컴퓨팅

10. Footnotes

1. 물론 공짜는 아니다. 원자를 고전 입자로 취급하는 순간 화학 결합의 생성·절단(전자 재배치)은 기본적으로 다루지 못한다. 반응을 보려면 반응성 힘장(ReaxFF)이나, 아예 매 스텝 밀도범함수이론으로 힘을 계산하는 제일원리 MD(ab initio MD)로 넘어가야 한다. 후자는 정확한 대신 계산비용이 살인적이다. ↩

2. 힘장 파라미터를 어떻게 정하느냐가 사실상 예술의 영역이다. 실험값과 양자계산을 짜맞춰 피팅하는데, 한 힘장에서 잘 맞던 파라미터가 다른 계로 가면 어긋나는 “전이성(transferability)” 문제가 MD판의 영원한 골칫거리다. ↩

3. 그래서 수소 관련 결합을 강체로 고정하는 SHAKE/LINCS 같은 구속 알고리즘을 써서 가장 빠른 진동을 제거하고, 시간 스텝을 2 fs로 늘리는 꼼수가 관행이다. LINCS 논문(Hess 1997)이 GROMACS 기본값의 출처다. 1 fs든 2 fs든, 어차피 우리가 보고 싶은 현상 앞에서는 티끌만 한 시간이다. ↩