본문으로 건너뛰기

DFT 이론 — 다체 문제에서 Kohn–Sham까지

이 튜토리얼의 모든 전자구조 계산 — SIESTA의 SCF(챕터 01)부터 VASP relax(챕터 03), TranSIESTA transport까지 — 은 밀도범함수론(density functional theory, DFT) 위에서 돌아간다. 이 장은 그 이론적 토대를 세운다. 다체 Schrödinger 방정식이 왜 직접 풀 수 없는 문제인지, Hohenberg–Kohn 정리가 무엇을 보장하는지, 그리고 Kohn–Sham 구성이 이 보장을 어떻게 실용적인 계산 절차로 바꾸는지를 정리한다. 이론 장이므로 입력 파일과 실행 섹션은 없다.

학습 목표

  • 다체 파동함수의 차원이 전자 수에 대해 지수적으로 커지는 이유와 Born–Oppenheimer 근사의 내용을 설명한다.
  • Hohenberg–Kohn 두 정리의 진술과 의미를 파악한다 — 밀도가 모든 것을 결정하고, 변분 원리가 성립한다.
  • Kohn–Sham 방정식과 유효 퍼텐셜 veff=vext+vH+vxcv_{\mathrm{eff}} = v_{\mathrm{ext}} + v_H + v_{xc}의 구조를 이해한다.
  • self-consistency 루프의 각 단계를 챕터 01 실습의 SCF 키워드(dDmax, mixing)와 대응시킨다.
  • XC functional 사다리(LDA → GGA → meta-GGA → hybrid)에서 각 단이 무엇을 추가하고 무엇을 대가로 치르는지 구분한다.
  • KS 고유값의 지위와 band gap 과소평가 문제의 기원을 개괄한다.

1. 다체 문제와 지수 폭발

원자핵과 전자로 이루어진 계의 비상대론적 Hamiltonian은 알려져 있다. 핵을 잠시 잊고 전자 NN개만 보면(원자 단위):

H^=12i=1Ni2+i=1Nvext(ri)+12ij1rirj\hat H = -\frac{1}{2}\sum_{i=1}^{N} \nabla_i^2 + \sum_{i=1}^{N} v_{\mathrm{ext}}(\mathbf r_i) + \frac{1}{2}\sum_{i \neq j} \frac{1}{|\mathbf r_i - \mathbf r_j|}

첫 항은 운동에너지, 둘째 항은 핵이 만드는 외부 퍼텐셜, 셋째 항이 전자–전자 Coulomb 반발이다. 방정식 자체는 완전하지만, 해인 다체 파동함수 Ψ(r1,,rN)\Psi(\mathbf r_1, \ldots, \mathbf r_N)3N3N개 변수의 함수라는 점이 문제다. 각 좌표축을 격자점 10개로만 표현해도 필요한 값의 개수는 103N10^{3N} — 원자 8개로 이루어진 탄소 사슬(전자 48개)이면 1014410^{144}개로, 저장 자체가 불가능하다.

지수 폭발을 피하는 명백한 시도 — Ψ\Psi를 1전자 함수들의 곱으로 쓰는 것 — 은 셋째 항 때문에 정확하지 않다. 전자–전자 상호작용이 좌표들을 얽어 놓아 변수 분리가 성립하지 않으며, 이 분리 불가능성이 상관(correlation)이다. 다체 문제를 정면으로 푸는 길은 막혀 있고, 우회로가 필요하다.

2. Born–Oppenheimer 근사

첫 번째 우회는 핵과 전자의 분리다. 핵은 전자보다 최소 세 자릿수 무겁다(양성자/전자 질량비 약 1836). 전자가 핵의 움직임에 즉각 적응한다고 보고, 핵 좌표를 고정한 채 전자 문제만 푼다 — 핵 위치 {R}\{\mathbf R\}는 방정식의 변수가 아니라 파라미터가 된다.

핵 배치마다 전자 바닥상태 에너지를 얻으면 그것이 퍼텐셜 에너지 곡면 E({R})E(\{\mathbf R\})이고, 그 기울기가 핵에 걸리는 힘이다. Hellmann–Feynman 정리 덕분에 이 힘은 파동함수를 핵 좌표로 미분하지 않고도 수렴된 전자 밀도로부터 직접 계산된다. 구조 최적화(챕터 03)가 "힘이 0인 핵 배치 찾기"로 정의되는 근거가 바로 이 근사다. 이후 이 장의 논의는 전부 고정된 핵 아래의 전자 문제다.

3. Hohenberg–Kohn 정리

전자 문제는 여전히 3N3N차원이다. Hohenberg와 Kohn(1964)은 여기서 급진적인 축약이 가능함을 증명했다. 두 정리를 증명 없이 진술한다.

정리 1 (밀도가 모든 것을 결정한다). 바닥상태 전자 밀도 n(r)n(\mathbf r)은 외부 퍼텐셜 vext(r)v_{\mathrm{ext}}(\mathbf r)을 (상수 차이를 제외하고) 유일하게 결정한다. 퍼텐셜이 정해지면 Hamiltonian이 정해지고 파동함수와 모든 관측량이 정해지므로, 바닥상태의 모든 성질은 원리적으로 밀도만의 범함수(functional)다. 3N3N개 변수의 Ψ\Psi 대신 3개 변수의 n(r)n(\mathbf r)이 기본 변수가 될 수 있다는 존재 보증이다.

정리 2 (변분 원리). 에너지 범함수 E[n]=F[n]+vext(r)n(r)drE[n] = F[n] + \int v_{\mathrm{ext}}(\mathbf r)\, n(\mathbf r)\, d\mathbf r가 존재하여, 임의의 허용 가능한 밀도에 대해 E[n]E0E[n] \geq E_0이고 등호는 참 바닥상태 밀도에서만 성립한다. 여기서 F[n]F[n]은 운동에너지와 전자–전자 상호작용을 담는 보편 범함수로, 어떤 계인지(vextv_{\mathrm{ext}}가 무엇인지)와 무관하게 같은 형태다.

의미를 분명히 해 두자. 두 정리는 "밀도만으로 충분하다"는 것을 보장할 뿐, F[n]F[n]의 구체적 형태를 알려주지 않는다 — 존재 증명이지 구성 방법이 아니다. F[n]F[n]을 정확히 아는 것은 원래의 다체 문제를 푸는 것과 같은 난이도다. DFT가 실용적 도구가 된 것은 다음 절의 Kohn–Sham 구성 덕분이다.

4. Kohn–Sham 구성 — 상호작용 없는 보조계

Kohn과 Sham(1965)의 아이디어는 이렇다. 참 계와 같은 바닥상태 밀도를 가지는, 상호작용 없는 가상의 보조계를 도입한다. 상호작용 없는 전자들의 운동에너지 Ts[n]T_s[n]과 고전적 Coulomb(Hartree) 에너지 EH[n]E_H[n]은 정확히 계산할 수 있으므로, 에너지 범함수를 다음과 같이 분해한다.

E[n]=Ts[n]+vextndr+EH[n]+Exc[n]E[n] = T_s[n] + \int v_{\mathrm{ext}}\, n\, d\mathbf r + E_H[n] + E_{xc}[n]

마지막 항 Exc[n]E_{xc}[n]교환–상관(exchange–correlation) 범함수다 — 참 운동에너지와 TsT_s의 차이, 그리고 고전적 Coulomb을 넘어서는 모든 상호작용 효과를 담는 "나머지"로 정의된다. 교환은 파동함수의 반대칭성(Pauli 원리)에서 오는 몫이고, 상관은 그 너머의 모든 것이다. 다체 문제의 모든 어려움이 이 한 항에 격리되었고, 알 수 없는 것은 이것뿐이다.

이 범함수에 정리 2의 변분 원리를 적용하면 1전자 형태의 방정식이 나온다.

[122+veff(r)]ψi(r)=εiψi(r),veff=vext+vH+vxc\left[-\frac{1}{2}\nabla^2 + v_{\mathrm{eff}}(\mathbf r)\right]\psi_i(\mathbf r) = \varepsilon_i\,\psi_i(\mathbf r), \qquad v_{\mathrm{eff}} = v_{\mathrm{ext}} + v_H + v_{xc}

이것이 Kohn–Sham 방정식이다. vH=δEH/δnv_H = \delta E_H/\delta n은 밀도가 만드는 고전적 Coulomb 퍼텐셜, vxc=δExc/δnv_{xc} = \delta E_{xc}/\delta n은 교환–상관 퍼텐셜이며, 밀도는 점유된 KS 궤도로부터 재구성된다.

n(r)=ioccψi(r)2n(\mathbf r) = \sum_i^{\mathrm{occ}} |\psi_i(\mathbf r)|^2

형태는 1전자 Schrödinger 방정식이지만, veffv_{\mathrm{eff}}가 밀도에 의존하므로 실질은 비선형 문제다. 그리고 ExcE_{xc}를 근사해야 하는 순간부터 DFT는 근사 이론이 된다 — 어떤 근사를 쓰는지가 6절의 주제다.

짚어 둘 것 하나. KS 궤도 ψi\psi_i와 고유값 εi\varepsilon_i는 형식적으로는 밀도를 만들기 위한 보조 객체다. 그럼에도 실무에서는 밴드 구조·DOS·PDOS(챕터 02)를 전부 KS 고유값으로 그리고, 대부분의 경우 실험 스펙트럼의 정성적 구조를 잘 재현한다. 다만 정량적 지위 — 특히 gap — 에는 구조적 한계가 있으며, 이것이 7절과 챕터 11의 주제다.

5. Self-consistency 루프

veffv_{\mathrm{eff}}nn에 의존하고, nnveffv_{\mathrm{eff}}의 해 ψi\psi_i로부터 나온다. 따라서 해는 고정점(fixed point)이고, 반복법으로 찾는다.

① 초기 밀도 n⁰ 추정 (보통 중성 원자 밀도의 중첩)
② v_eff[n] 구성: Hartree + XC 퍼텐셜 계산
③ KS Hamiltonian 대각화 → 궤도 {ψ_i}와 고유값 {ε_i}
④ 점유 궤도로 새 밀도 n_out 재구성
⑤ 수렴 판정: n_out과 n_in의 차이가 tolerance 이하인가?
예 → 종료. 총에너지·힘·고유값 산출
아니오 → ⑥ mixing: n_in(new) = n_in과 n_out의 혼합 → ②로 복귀

Kohn–Sham SCF 루프 순서도

그림 1. self-consistency 루프(개념도) — 밀도로 H[ρ]H[\rho]를 만들고, KS 방정식을 풀어 새 밀도를 얻고, dDmax가 tolerance 아래로 내려올 때까지 mixing으로 반복한다.

⑥에서 noutn_{\mathrm{out}}을 통째로 다음 입력으로 쓰지 않고 일부만 섞는 이유는 안정성이다. 금속처럼 유전 응답이 긴 파장에서 발산하는 계는 밀도의 작은 변화가 퍼텐셜의 큰 변화로 증폭되어(charge sloshing) 단순 반복이 진동·발산한다. 새 밀도를 조금씩 반영하고(선형 mixing) 이전 반복들의 이력으로 최적 혼합을 추정하는(Pulay mixing) 것이 표준 처방이다.

이 루프의 각 단계는 챕터 01 실습의 키워드·출력과 1:1로 대응한다.

루프 단계SIESTA 키워드 / 출력 (챕터 01)
① 초기 밀도DM.UseSaveDM — 저장된 밀도 행렬이 있으면 그것으로 시작 (warm start)
③ 대각화SCF 1회 반복의 주 비용. iscf 카운터가 반복 횟수
⑤ 수렴 판정dDmax (밀도 행렬 원소의 반복 간 최대 변화)가 SCF.DM.Tolerance 이하
⑥ mixingSCF.Mixer.Weight 0.2 (새 밀도 20% 반영), SCF.Mixer.History 15 (Pulay 이력)

수렴 실패는 두 형태로 나타난다. dDmax가 줄지 않고 진동하면 반복 지도가 불안정한 것이므로 mixing weight를 줄이는 것이 첫 번째 처방이다 — 반복의 보폭을 줄여 안정 영역 안에 머무르게 하는 것이다. 반대로 MaxSCFIterations에 걸려 멈춘 계산은 수렴한 것이 아니며 결과를 쓰면 안 된다는 판정 규칙(챕터 01)도 이 구조에서 나온다.

Eharris와 E_KS — 챕터 01 출력의 두 에너지 열

챕터 01 stdout의 EharrisE_KS는 같은 총에너지에 대한 두 추정치다. Harris 범함수(Eharris)는 그 반복의 입력 밀도만으로 평가되고, Kohn–Sham 범함수(E_KS)는 출력 밀도로 평가된다. SCF가 수렴해 입력·출력 밀도가 일치하면 두 값은 하나로 합쳐진다. 따라서 dDmax와 함께 두 열의 차이가 줄어드는 것 자체를 수렴 지표로 읽을 수 있다.

6. XC functional 사다리

남은 것은 Exc[n]E_{xc}[n]의 근사다. 근사들은 밀도로부터 얼마나 많은 정보를 쓰는지에 따라 사다리처럼 배열된다.

대표입력 정보특징과 대가
LDAPZ, PW92국소 밀도 n(r)n(\mathbf r)균일 전자 기체의 ϵxc(n)\epsilon_{xc}(n)을 각 지점에 적용. 놀랍도록 견고하지만 결합을 과대평가(격자상수 과소)하는 경향
GGAPBEnn, n\nabla n밀도 기울기 보정. 결합 에너지·구조 크게 개선, 비용은 LDA와 사실상 동일. 고체 계산의 사실상 표준
meta-GGASCAN 등nn, n\nabla n, 운동에너지 밀도 τ\tau결합 유형 구분 능력 향상. 비용 소폭 증가, 수치 안정성에 주의
hybridB3LYP, HSE06, PBE0+ 비국소 HF 교환 일부(보통 20–25%) 혼합self-interaction 오차 완화, gap·준위 위치 개선. 비국소 교환 때문에 비용이 한두 자릿수 증가

이 위계를 Perdew는 "Jacob의 사다리"라 불렀다 — 위로 갈수록 대체로 정확해지지만 비용이 커지고, 어떤 단도 모든 문제에서 아래 단을 이기지는 못한다. 실무 규칙 하나가 여기서 나온다: 비교하려는 양들은 반드시 같은 functional로 계산한다. 서로 다른 단에서 얻은 에너지를 직접 빼는 것은 의미가 없다.

이 튜토리얼은 PBE를 표준으로 쓴다. 이유는 두 가지다. 첫째, SIESTA와 VASP 양쪽에서 동일하게 지원되므로 LCAO/plane-wave 교차 검증(챕터 03)과 transport 워크플로우 전체를 같은 functional로 관통할 수 있다. 둘째, 구조·총에너지·힘에 대해서는 PBE가 비용 대비 충분히 정확하다.

hybrid는 전체 transport 계산에 쓰기에는 비용이 지나치게 크고 LCAO NEGF 구현과의 결합도 제한적이므로, 고립 분자의 준위 기준점(reference) 산출에만 쓴다 — 그 용법이 챕터 11의 DFT+Σ 보정이다.

7. Band gap 과소평가 — 예고

한 가지 한계는 미리 알아 두어야 한다. KS 고유값 εi\varepsilon_i는 보조계의 라그랑주 승수이지 전자를 실제로 넣고 뺄 때의 에너지(준입자 에너지)가 아니다. semi-local functional(LDA/GGA)은 self-interaction 오차와 derivative discontinuity 누락 때문에 점유 준위를 너무 높게, 비점유 준위를 너무 낮게 놓는다 — 그 결과 band gap이 체계적으로 과소평가되며, 전형적 반도체에서 수십 % 오차가 흔하다.

밴드 개형이나 구조 예측에는 큰 문제가 아니지만, transport에서는 치명적일 수 있다. 공명 준위의 위치가 수백 meV 틀리면 gap 안의 T(EF)T(E_F)는 자릿수 단위로 왜곡되기 때문이다. 이 오차의 기원과 문헌 표준 처방(DFT+Σ)은 챕터 11에서 다룬다.

핵심 요약

개념내용
지수 폭발Ψ(r1,,rN)\Psi(\mathbf r_1,\ldots,\mathbf r_N)3N3N차원 — 직접 풀이 불가
Born–Oppenheimer핵 고정, 전자 문제만. E({R})E(\{\mathbf R\})가 퍼텐셜 곡면, 그 기울기가 힘
HK 정리 1바닥상태 밀도가 vextv_{\mathrm{ext}}와 모든 성질을 결정
HK 정리 2E[n]E0E[n] \geq E_0, 등호는 참 밀도 — 변분 원리
KS 구성같은 밀도의 비상호작용 보조계. veff=vext+vH+vxcv_{\mathrm{eff}} = v_{\mathrm{ext}} + v_H + v_{xc}
SCF 루프밀도 → 퍼텐셜 → 대각화 → 새 밀도 → mixing → 수렴 (dDmax, SCF.Mixer.*)
Eharris vs E_KS입력 밀도 / 출력 밀도로 평가한 총에너지 — 수렴 시 일치
XC 사다리LDA → GGA(PBE) → meta-GGA → hybrid. 이 튜토리얼은 PBE + hybrid reference
KS gap준입자 gap이 아님 — semi-local에서 체계적 과소평가 (챕터 11)

연습문제

  1. 지수 폭발의 감각. 좌표축당 격자점 10개로 파동함수를 표로 저장한다고 하자. 전자 2개(He 원자)와 전자 48개(C 사슬 8원자)에 필요한 저장 용량을 각각 추정하고, 밀도 n(r)n(\mathbf r)을 같은 격자에 저장할 때와 비교하라.
  2. Hartree 퍼텐셜. EH[n]=12n(r)n(r)rrdrdrE_H[n] = \frac{1}{2}\iint \frac{n(\mathbf r)\,n(\mathbf r')}{|\mathbf r - \mathbf r'|}\, d\mathbf r\, d\mathbf r'의 범함수 미분이 vH(r)=n(r)rrdrv_H(\mathbf r) = \int \frac{n(\mathbf r')}{|\mathbf r - \mathbf r'|}\, d\mathbf r'임을 보여라.
  3. mixing의 필요성. SCF 반복을 고정점 방정식 n=g(n)n = g(n)의 반복법으로 보면, 수렴 조건은 반복 지도의 기울기 크기가 1보다 작은 것이다. 선형 mixing n(k+1)=(1α)n(k)+αg(n(k))n^{(k+1)} = (1-\alpha)\,n^{(k)} + \alpha\, g(n^{(k)})이 왜 발산하는 반복을 수렴시킬 수 있는지, 그리고 α\alpha를 무작정 줄이면 무엇을 대가로 치르는지 설명하라.
  4. Harris vs KS 에너지. 챕터 01 계산의 stdout에서 첫 반복과 마지막 반복의 EharrisE_KS 차이를 비교하라. 왜 반복이 진행될수록 두 값이 접근하는지 5절의 정의로 설명하라.
  5. KS 고유값의 의미. 6절과 7절을 바탕으로, "PBE 밴드 구조의 gap이 실험보다 작다"는 관찰이 계산 실수가 아니라 이론의 구조적 한계인 이유를 두 문장으로 정리하라. KS 방정식을 실제로 푸는 수치 골격 — 기저 전개와 k-샘플링 — 은 다음 장 — Bloch 정리와 기저에서 이어진다.

Ref: Hohenberg & Kohn, Phys. Rev. 136, B864 (1964); Kohn & Sham, Phys. Rev. 140, A1133 (1965); R. M. Martin, Electronic Structure: Basic Theory and Practical Methods (Cambridge University Press, 2004).