움직이는 교과서 · Interactive Textbook

Python & DFT 기초

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

Python ASE GPAW DFT Bulk Properties Convergence

Section 1.1

계산재료과학 소개

기후변화에 대응하기 위해 효율적인 촉매 개발의 중요성이 커지고 있습니다. 밀도범함수이론(Density Functional Theory, DFT)은 원자핵 주변의 전자 밀도를 시뮬레이션하여 촉매의 활성점, 흡착 상호작용, 전자구조 등의 특성을 예측할 수 있게 해줍니다.

계산재료과학에서는 다양한 시뮬레이션 방법이 사용되며, 각 방법은 고유한 길이 스케일과 시간 스케일 영역을 다룹니다. DFT는 원자 수준(~1 nm)에서 전자구조를 직접 계산하고, 분자동역학(MD)은 더 큰 시스템의 원자 운동을, kinetic Monte Carlo(kMC)와 연속체 모델은 더 거시적인 현상을 다룹니다.

$$H\Psi = E\Psi \quad \xrightarrow{\text{DFT}} \quad E[n(\mathbf{r})] = T_s[n] + V_{\text{ext}}[n] + E_H[n] + E_{xc}[n]$$
Figure 1.1 Multi-Scale Simulation Map
관심 시스템 크기 1 nm
적합한 방법 DFT
원자 수 범위 ~100
시간 스케일 ~ps
DFT는 원자 수준의 전자구조를 다루는 가장 기본적인 방법입니다
시스템이 커지면 DFT로는 계산 비용이 너무 높아지므로, 더 큰 스케일에서는 MD, kMC, 연속체 모델 등으로 넘어가야 합니다. 각 방법의 적용 범위를 이해하는 것이 멀티스케일 시뮬레이션의 핵심입니다.
  1. 슬라이더를 왼쪽 끝(~0.1 nm)으로 이동하세요. DFT 영역이 강조됩니다.
  2. 천천히 오른쪽으로 이동하면서 어떤 시뮬레이션 방법이 적합한지 관찰하세요.
  3. ~10 nm 근처에서 DFT에서 MD 영역으로 전환되는 것을 확인하세요.
  4. ~1 μm 이상이 되면 kMC나 연속체 모델이 필요해집니다.
Key Insight: 시뮬레이션 방법의 선택은 연구 대상의 길이/시간 스케일에 의해 결정됩니다. DFT는 가장 정밀하지만 가장 비싸고, 연속체 모델은 가장 빠르지만 원자 수준의 정보를 잃습니다.

Section 1.2

Python & ASE 기초

ASE(Atomic Simulation Environment)는 원자 시뮬레이션을 위한 Python 라이브러리입니다. ASE의 핵심은 Atoms 객체와 Calculator 패턴입니다. Atoms 객체에 원자의 종류, 위치, 셀 정보를 담고, Calculator(예: GPAW)를 연결하면 에너지, 힘 등을 계산할 수 있습니다.

가장 간단한 예로 H2 분자의 포텐셜을 Morse potential로 모델링할 수 있습니다. Morse potential은 결합 에너지($D_e$), 평형 결합 거리($r_e$), 포텐셜 폭 파라미터($\alpha$)로 완전히 결정됩니다.

$$V(r) = D_e\bigl(1 - e^{-\alpha(r - r_e)}\bigr)^2 - D_e$$
Figure 1.2 Morse Potential for H2
목표: De를 최대로, re를 최소로 설정했을 때 물리적으로 어떤 결합을 의미할까요?
짧은 결합 거리에 큰 해리 에너지 — 이것은 매우 강한 공유결합을 의미합니다. 예를 들어 N2 삼중결합이 이런 특성을 보입니다. 슬라이더를 움직여 readout의 zero-point frequency가 어떻게 변하는지도 관찰해보세요.
현재: 슬라이더를 조작해보세요
De 4.75 eV
re 0.74 Å
α 1.94 Ź
Equilibrium r 0.74 Å
Dissociation E 4.75 eV
ZP Frequency —
V(r) potential
F(r) = −dV/dr
$r = r_e$에서 힘이 0이 되고, 포텐셜이 최소입니다
$r < r_e$에서는 척력(양의 힘)이, $r > r_e$에서는 인력(음의 힘)이 작용합니다. 이 균형점이 평형 결합 거리이며, ASE에서 구조 최적화는 바로 이 점을 찾는 과정입니다.
  1. De = 4.75 eV, re = 0.74 Å, α = 1.94 Å−1로 시작하세요. 이것은 H2의 실험 파라미터입니다.
  2. ASE에서 Atoms('H2', positions=[(0,0,0),(0,0,d)])로 H2를 정의합니다. d가 r에 해당합니다.
  3. Calculator를 연결하면 atoms.get_forces()로 각 원자에 작용하는 힘을 얻습니다.
  4. α를 변경하면 포텐셜 우물의 폭이 변합니다. α가 클수록 우물이 좁아지고 진동 주파수가 높아집니다.
  5. De를 변경하면 우물의 깊이가 변합니다. 깊을수록 결합이 강합니다.
Key Insight: ASE의 Atoms + Calculator 패턴은 어떤 포텐셜(Morse, EMT, DFT 등)이든 동일한 인터페이스로 에너지와 힘을 계산합니다. 포텐셜만 바꾸면 같은 코드로 다른 수준의 계산이 가능합니다.

Section 1.3

K-point 수렴 테스트

주기적 시스템에서 Bloch 정리에 의해 전자 상태는 Brillouin zone 내의 k-point들로 레이블됩니다. 실제 계산에서는 유한한 수의 k-point로 Brillouin zone을 샘플링해야 합니다. k-point가 너무 적으면 에너지에 큰 오차가 발생하고, 너무 많으면 불필요한 계산 비용이 들어갑니다.

아래 그래프는 Al(FCC) bulk에 대한 k-point 수렴 테스트를 시뮬레이션한 것입니다. k-point 수를 늘리면 에너지가 특정 값으로 수렴하는 것을 볼 수 있습니다. 수렴 기준(convergence threshold)을 설정하면 필요한 최소 k-point 수를 결정할 수 있습니다.

$$E(k) = E_{\text{conv}} + A\,e^{-Bk} + C\sin\!\left(\frac{\pi k}{2}\right)e^{-k/3}$$
Figure 1.3 K-point Convergence (Al FCC)
수렴 기준 0.010 eV
수렴 k-point —
에너지 차이 —
수렴 에너지 −4.040 eV
k-point가 충분히 크면 에너지가 수렴합니다
금속 시스템(예: Al)은 Fermi level 근처의 전자 상태가 k-space에서 빠르게 변하므로 많은 k-point가 필요합니다. 반도체나 절연체는 상대적으로 적은 k-point로도 수렴합니다.
Q: k-point를 무한히 늘리면 계산 정확도는 어떻게 되나요?
k-point를 늘리면 Brillouin zone 적분의 정밀도가 높아지지만, 충분히 많아지면 추가적인 개선이 무시할 수 있을 정도로 작아져 에너지가 수렴합니다. 이것이 수렴 테스트가 필요한 이유입니다.

Section 1.4

Plane-wave Cutoff 수렴

Plane-wave DFT에서 Kohn-Sham 궤도함수(orbital)는 평면파(plane wave)의 합으로 전개됩니다. Cutoff energy($E_{\text{cut}}$)는 이 전개에 포함되는 평면파의 최대 운동에너지를 결정합니다. 높은 cutoff은 더 많은 평면파를 포함하여 급격한 전자 밀도 변화를 잘 표현하지만, 계산 비용도 증가합니다.

GPAW에서는 mode=PW(cutoff)으로 cutoff energy를 설정합니다. 아래 그래프는 cutoff에 따른 에너지 수렴을 보여줍니다. 일반적으로 cutoff이 증가하면 에너지가 단조감소하며 수렴합니다.

$$\psi_{n\mathbf{k}}(\mathbf{r}) = \sum_{|\mathbf{G}+\mathbf{k}|^2 \le 2E_{\text{cut}}} c_{n\mathbf{k}+\mathbf{G}}\,e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}}$$
Figure 1.4 Cutoff Energy Convergence (Al FCC)
Econv −4.04 eV
A (amplitude) 1000
수렴 Cutoff —
에너지 차이 —
상대 계산 비용 —
Cutoff이 높을수록 정확하지만 계산 비용이 급격히 증가합니다
Plane-wave 수는 $E_{\text{cut}}^{3/2}$에 비례하므로, cutoff을 2배로 올리면 basis set 크기가 약 2.8배 증가합니다. 따라서 수렴한 최소 cutoff을 사용하는 것이 효율적입니다.
  1. Econv = −4.04 eV, A = 1000으로 시작하세요.
  2. 그래프에서 cutoff 100 eV의 에너지와 800 eV의 에너지 차이를 비교하세요.
  3. A를 2000으로 올려보세요. 낮은 cutoff에서의 오차가 더 커집니다.
  4. A를 500으로 줄여보세요. 낮은 cutoff에서도 수렴이 빨라집니다. 이것은 원소에 따라 다릅니다 — 산소 같은 1주기 원소는 급격한 전자 밀도를 가져 높은 cutoff이 필요합니다.
Key Insight: GPAW에서 mode=PW(300)은 300 eV cutoff을 의미합니다. 실제 연구에서는 원소와 pseudopotential에 따라 400–600 eV 사이의 cutoff이 흔히 사용됩니다.

Section 1.5

상태방정식 (Equation of State)

DFT 계산에서 격자 상수를 결정하기 위해 Birch-Murnaghan 상태방정식(EOS)을 사용합니다. 여러 격자 상수(또는 부피)에서 에너지를 계산한 뒤, E(V) 곡선을 Birch-Murnaghan 식으로 fitting하면 평형 격자 상수, 최소 에너지, 체적탄성률(bulk modulus)을 정밀하게 추출할 수 있습니다.

FCC 구조에서 부피 $V = a^3/4$ (primitive cell)이므로, 격자 상수 $a$를 변화시키면 부피가 변합니다. 아래 그래프에서 슬라이더를 조절하여 Al(FCC)의 EOS를 탐색해보세요.

$$E(V) = E_0 + \frac{9V_0 B_0}{16}\left\{\left[\left(\frac{V_0}{V}\right)^{2/3}-1\right]^3 B_0' + \left[\left(\frac{V_0}{V}\right)^{2/3}-1\right]^2\left[6-4\left(\frac{V_0}{V}\right)^{2/3}\right]\right\}$$
유도 과정 보기
1
Euler strain: $f = \frac{1}{2}\left[\left(\frac{V_0}{V}\right)^{2/3} - 1\right]$
유한 변형(finite strain)을 정의합니다. $V = V_0$일 때 $f = 0$입니다.
2
에너지를 strain의 Taylor 급수로 전개: $E(f) = E_0 + \alpha f^2 + \beta f^3 + \cdots$
평형 근처에서 $f$의 다항식으로 에너지를 근사합니다.
3
$B_0 = -V\frac{dP}{dV}\big|_{V_0}$, $B_0' = \frac{dB}{dP}\big|_{P=0}$
Bulk modulus $B_0$는 압축에 대한 저항, $B_0'$는 압력에 따른 bulk modulus의 변화율입니다.
4
3차까지의 계수를 $B_0$, $B_0'$로 표현하면 Birch-Murnaghan EOS가 됩니다.
3개의 피팅 파라미터($E_0$, $V_0$, $B_0$, $B_0'$)로 E-V 곡선을 정확히 기술합니다.
Figure 1.5 Birch-Murnaghan Equation of State (Al FCC)
목표: 실험값 a = 4.05 Å에 가장 가까운 B0, B0'를 찾아보세요
Al의 실험 bulk modulus는 약 76 GPa, B0'는 약 4.2입니다. 슬라이더를 조절해서 실험값에 가까운 EOS 곡선을 만들어보세요. 곡선의 모양이 어떻게 변하는지 관찰하세요.
현재: 슬라이더를 조작해보세요
a0 4.05 Å
B0 76 GPa
B0' 4.2
Optimal a —
Bulk Modulus —
Min Energy —
EOS fitting으로 격자 상수를 정밀하게 결정합니다
단순히 몇 개의 격자 상수에서 에너지를 비교하는 것보다, Birch-Murnaghan fitting을 사용하면 연속적인 E(V) 곡선에서 최소점을 찾으므로 더 정확한 평형 격자 상수를 얻을 수 있습니다.
Q: Birch-Murnaghan EOS에서 B0는 무엇을 나타내나요?
B0는 체적탄성률(bulk modulus)로, 재료가 균일한 압축에 저항하는 정도를 나타냅니다. 단위는 GPa이며, B0 = −V(dP/dV)로 정의됩니다. 값이 클수록 재료가 단단합니다.

What If?

만약에…

만약 k-point를 1로 줄이면?
Brillouin zone의 Γ점만 샘플링하게 되어 금속에서는 큰 오차가 발생합니다. 특히 Fermi level 근처의 전자 상태를 제대로 반영하지 못해 에너지, 힘, 전자구조 모두 부정확해집니다.
클릭해서 확인 →
만약 cutoff을 50 eV로 낮추면?
Basis set이 너무 작아 전자의 급격한 밀도 변화를 표현하지 못합니다. 결과적으로 에너지가 크게 틀어지고, 원자간 힘도 부정확해져 구조 최적화가 잘못된 결과를 줍니다.
클릭해서 확인 →
만약 격자상수를 실험값 대신 임의로 넣으면?
잘못된 평형 구조에서 출발하면 모든 후속 계산(표면에너지, 흡착에너지 등)이 체계적으로 틀어집니다. EOS fitting으로 정확한 격자 상수를 결정한 후 본 계산을 시작해야 합니다.
클릭해서 확인 →

Key Takeaways

Chapter 1 핵심 메시지

✓
수렴 테스트 필수
k-point와 cutoff은 반드시 수렴을 확인한 후 본 계산에 사용해야 합니다. 수렴하지 않은 파라미터로 계산하면 모든 결과가 부정확해집니다.
EOS
EOS로 격자 상수 결정
Birch-Murnaghan fitting으로 정확한 평형 격자 상수와 체적탄성률을 얻을 수 있습니다. 이것이 모든 후속 표면/흡착 계산의 출발점입니다.
⇌
정확도 vs 비용
더 높은 k-point와 cutoff은 정확하지만 계산 비용이 급격히 증가합니다. 수렴한 최소 파라미터를 찾아 효율적으로 계산하는 것이 핵심입니다.