움직이는 교과서 · Interactive Textbook

분자동역학 시뮬레이션

Choung, Simulation Tutorials for Computational Materials Science, Chapter 3. 슬라이더를 움직이면 그래프가 실시간으로 반응합니다.

Molecular Dynamics Lennard-Jones NPT Thermal Expansion Diffusion

Section 3.1

분자동역학의 원리

분자동역학(Molecular Dynamics, MD) 시뮬레이션은 원자 사이의 상호작용 퍼텐셜로부터 힘을 계산하고, Newton의 운동 방정식 $F = ma = -dV/dr$을 수치적으로 적분하여 원자의 궤적을 추적하는 방법입니다.

시간 적분에는 Verlet algorithm이 널리 사용됩니다. 이 방법은 시간 가역성(time-reversibility)을 보장하고, 에너지 보존 성능이 뛰어나며, 구현이 간단합니다. 원자 간 상호작용을 기술하는 가장 기본적인 모델이 Lennard-Jones potential입니다.

$$F = ma = -\frac{dV}{dr}$$
유도 과정 보기: Verlet Integration
1
Taylor 전개: $r(t+\Delta t) = r(t) + v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 + \frac{1}{6}\dot{a}(t)\Delta t^3 + \cdots$
현재 위치에서 미래 위치를 급수로 전개합니다.
2
반대 방향으로도 전개: $r(t-\Delta t) = r(t) - v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 - \frac{1}{6}\dot{a}(t)\Delta t^3 + \cdots$
과거 위치도 같은 방식으로 구합니다.
3
두 식을 더하면 홀수 차수 항이 상쇄: $r(t+\Delta t) + r(t-\Delta t) = 2r(t) + a(t)\Delta t^2 + O(\Delta t^4)$
속도 항이 자연스럽게 소거되어, 속도 없이도 위치를 업데이트할 수 있습니다.
4
최종 Verlet formula: $r(t+\Delta t) = 2r(t) - r(t-\Delta t) + a(t)\Delta t^2$
가속도 $a(t) = F(t)/m$은 퍼텐셜의 기울기로부터 계산합니다. 오차는 $O(\Delta t^4)$으로, 매우 정확합니다.
Figure 3.1 Lennard-Jones Potential
목표: rmin에서 force가 정확히 0이 되는지 확인해보세요
빨간색 force 곡선이 0을 지나는 위치와, 검은색 potential 곡선의 최솟값 위치가 일치하는지 확인해보세요. 이것은 $F = -dV/dr$의 직접적인 결과입니다.
현재: 슬라이더를 조작해보세요
ε (well depth) 0.167 eV
σ 2.55 Å
rmin --
Well Depth --
F=0 Distance --
V(r) (eV)
F(r) (eV/Å)
  1. ε을 0.167 eV (Cu EMT 기본값)으로 설정하고, σ를 2.55 Å로 맞추세요.
  2. rmin = 21/6σ ≈ 2.86 Å에서 potential이 최솟값(-ε)을 가지는 것을 확인하세요.
  3. 같은 위치에서 force(빨간색)가 정확히 0을 지나는 것을 확인하세요. 이것이 평형 거리입니다.
  4. r < rmin: force가 양수(repulsive), r > rmin: force가 음수(attractive)임을 확인하세요.
Key Insight: Verlet algorithm은 이 force를 매 time step마다 계산하여 원자의 가속도를 구하고, 위치를 업데이트합니다. LJ potential의 기울기가 급한 영역(r < σ)에서는 작은 time step이 필수적입니다.

Section 3.2

Maxwell-Boltzmann 속도 분포

MD 시뮬레이션에서 원자의 초기 속도는 Maxwell-Boltzmann 분포에서 무작위로 추출합니다. 이 분포는 통계역학에서 유도되며, 주어진 온도에서 원자가 가질 수 있는 속력의 확률 분포를 나타냅니다.

온도가 높을수록 분포가 넓어지고 가장 확률이 높은 속력(most probable speed)이 커집니다. 반대로 질량이 큰 원자는 같은 온도에서도 느리게 움직입니다.

$$f(v) = 4\pi \left(\frac{m}{2\pi k_BT}\right)^{3/2} v^2 \exp\!\left(-\frac{mv^2}{2k_BT}\right)$$
Figure 3.2 Maxwell-Boltzmann Distribution
Temperature 300 K
Mass 63.5 amu
vp (most probable) --
Mean Speed --
vrms --
Key Observation
온도가 높을수록 분포가 넓어지고, 질량이 클수록 좁아집니다. Cu(63.5 amu)와 H(1 amu)를 비교해보면 속도 분포의 차이가 극적입니다.
  1. T = 300 K, m = 63.5 amu (Cu)에서 vp를 확인하세요.
  2. T를 600 K로 올려보세요. vp가 약 √2배 커지는 것을 확인할 수 있습니다.
  3. m을 1 amu (H)로 바꿔보세요. 같은 온도에서 수소는 구리보다 약 8배 빠릅니다.
  4. MD에서 time step은 가장 빠른 원자의 진동 주기에 맞춰야 합니다. H가 포함된 시스템에서 time step이 더 작아야 하는 이유를 이해할 수 있습니다.
Key Insight: Maxwell-Boltzmann 분포는 $v^2 \exp(-mv^2/2k_BT)$의 곱으로, 저속 영역은 위상 공간($v^2$)이, 고속 영역은 Boltzmann factor가 지배합니다.

Section 3.3

NPT 앙상블과 열평형

MD 시뮬레이션에서 앙상블(ensemble)은 어떤 열역학적 변수를 고정할지를 결정합니다. NVE(microcanonical)는 에너지가 보존되고, NVT(canonical)는 thermostat으로 온도를 제어하며, NPT(isothermal-isobaric)는 온도와 압력 모두를 제어합니다.

실험 조건(일정한 온도, 대기압)을 모사하려면 NPT 앙상블이 적합합니다. 열팽창이나 상전이 같은 열역학적 물성을 계산하려면 부피가 자유롭게 변할 수 있는 NPT가 필수적입니다.

$$T(t) = T_\text{target} + (T_\text{init} - T_\text{target})\exp(-t/\tau) + \text{noise}$$
Figure 3.3 MD Temperature Evolution
Ttarget 300 K
Tinit 10 K
τ (coupling) 1.5 ps
N (atoms) 64
T (at 10 ps) --
Equilibration Time --
Fluctuation ΔT --
Key Observation
Thermostat coupling time τ가 작을수록 빠르게 수렴하지만, 원자 수 N이 적으면 온도 요동이 커집니다. 실제 MD에서는 충분한 시스템 크기와 적절한 τ의 균형이 중요합니다.

Section 3.4

열팽창 계수

NPT MD 시뮬레이션에서 여러 온도에서의 격자 상수(lattice parameter)를 측정하면 열팽창 계수 $\alpha$를 구할 수 있습니다. 선형 열팽창 모델에서 $a(T) = a_0[1 + \alpha(T - T_\text{ref})]$이며, 실제로는 고온에서 비선형 기여($\beta$ 항)가 있습니다.

Notebook 05에서 Cu의 EMT potential을 사용한 NPT MD 결과를 NIST 실험 데이터와 비교한 바 있습니다. 시뮬레이션과 실험의 열팽창 계수를 비교하여 force field의 정확도를 평가할 수 있습니다.

$$a(T) = a_0\left[1 + \alpha(T - T_\text{ref}) + \beta(T - T_\text{ref})^2\right]$$
Figure 3.4 Thermal Expansion
목표: α를 조절하여 실험 데이터와 가장 잘 맞는 값을 찾아보세요
빨간색 점선(실험)과 검은색 실선(시뮬레이션)이 겹치도록 α 슬라이더를 조절하세요. 오차가 5% 이내이면 성공입니다.
현재: 슬라이더를 조작해보세요
a0 3.615 Å
α (×10-5/K) 1.70
αsim --
αexp --
Error --
Simulation
Experiment (NIST Cu)

Section 3.5

확산 계수와 MSD

Mean Square Displacement(MSD)는 원자가 시간에 따라 얼마나 멀리 이동했는지를 정량화합니다. Einstein relation에 의하면 확산 계수 $D$는 MSD의 장시간 기울기로부터 구할 수 있습니다: $\text{MSD} = 6Dt$ (3차원).

짧은 시간에서는 원자가 아직 주변 원자와 충돌하기 전이므로 탄도적(ballistic) 운동을 하며, MSD $\propto t^2$입니다. 충분한 시간이 지나면 확산(diffusive) regime으로 전환되어 MSD $\propto t$가 됩니다.

$$D = \lim_{t\to\infty}\frac{\text{MSD}(t)}{6t} = \lim_{t\to\infty}\frac{\langle|r(t)-r(0)|^2\rangle}{6t}$$
Figure 3.5 MSD Analysis
D (×10-9 m2/s) 2.3
Temperature 500 K
D (from fit) --
Crossover Time --
MSD (total)
Linear fit (diffusive)
Ballistic regime
Key Observation
확산 계수는 MSD의 기울기에서 직접 추출할 수 있습니다. D와 T를 조절하여 탄도-확산 전환 지점이 어떻게 변하는지 관찰해보세요.
만약 time step을 10 fs로 크게 잡으면?
빠르게 진동하는 경량 원자(H)의 운동을 놓쳐 에너지가 보존되지 않고 시뮬레이션이 불안정해집니다. 가장 빠른 진동의 주기 대비 time step이 너무 크면 numerical integration error가 축적됩니다.
클릭해서 확인 →
만약 NPT 대신 NVE로 돌리면?
온도와 압력이 제어되지 않아 열팽창 측정이 불가능합니다. NVE에서는 총 에너지만 보존되므로, 부피 변화가 허용되지 않고 온도도 drift합니다.
클릭해서 확인 →
만약 시뮬레이션 시간이 1 ps로 짧으면?
열평형에 도달하기 전이라 측정된 물성이 신뢰할 수 없습니다. 통계적 수렴을 위해 equilibration 시간 이후에도 충분한 production run이 필요합니다.
클릭해서 확인 →
Q: MSD가 시간에 선형적이면 운동 regime은?
MSD ∝ t (선형)이면 확산 regime입니다. 탄도적 regime에서는 MSD ∝ t²이고, 구속 상태에서는 MSD가 상수로 수렴합니다. Einstein relation D = MSD/(6t)는 확산 regime에서만 유효합니다.
Q: Maxwell-Boltzmann 분포의 most probable speed는 온도가 2배가 되면?
$v_p = \sqrt{2k_BT/m}$이므로, T가 2배이면 $v_p$는 $\sqrt{2}$배가 됩니다. 운동 에너지 $\frac{1}{2}mv^2 = \frac{3}{2}k_BT$에서 속력은 온도의 제곱근에 비례합니다.

Key Takeaways

Chapter 3 핵심 메시지

LJ
Lennard-Jones → 실제
LJ potential은 MD의 기본이지만, 실제 촉매 연구에서는 DFT나 MLP(Machine-Learned Potential)를 사용합니다. EMT도 금속에 한정된 간단한 모델입니다.
NPT
앙상블 선택
NPT는 실험 조건(일정 T, P)을 모사하므로 열팽창 등 열역학적 물성 계산에 적합합니다. NVE는 검증용, NVT는 표면/나노입자에 주로 사용합니다.
∞
충분한 시뮬레이션
열평형 + 통계적 수렴을 위해 충분한 시간과 시스템 크기가 필요합니다. 짧은 시뮬레이션은 비평형 상태의 물성을 측정하는 오류를 범합니다.