1. 개요
Density Functional Theory(DFT)[math(N)]개의 전자로 이루어진 분자의 파동함수를 표현하는데는 [math(4N)]개의 좌표가 필요하다.[1] 하지만 [math(N)]이 커질수록 파동함수는 우리가 필요한 정보 이상을 가지기도 하고 직관적으로 물리적 의미를 파악하기도 어려워진다. 그렇기에 파동함수 보다 더 적은 변수를 가지는 함수를 이용해 에너지나 다른 변수들의 기댓값을 구하려는 시도가 계속되었고, 그중 가장 유명한 이론이 전자 밀도를 이용한 밀도범함수 이론(Density Functional Theory, DFT)이다. 후술하겠지만 바닥 상태의 전자 밀도를 정확하게 안다면 신기하게도 계의 모양[2]이 하나로 정해지며 원론적으론 파동함수를 계산하지 않고도 모든 물리량의 기댓값을 계산할 수 있다! 물론 원론적으로만 그렇고 실제 전자 밀도를 실험 없이 바로 얻어낼 방법은 없으므로 이를 근사적으로 얻을 방법을 찾는 것이 DFT의 핵심이다.
2. 전자 확률 밀도
어떻게 파동함수와 전자 확률 밀도를 연관지을 수 있을까? 우리가 찾고자 하는 양은 점 [math(\textbf{r})]에서 미소부피 [math(d\textbf{r}(=dxdydz))] 안에서 전자를 발견할 확률이다. 그리고 1번 전자가 [math(m_{s1})] 스핀을 가지고 좌표 [math(\textbf{r}_1)]에서 미소부피 [math(d\textbf{r}_1)]안에서 발견될 확률, 2번 전자가 [math(m_{s2})] 스핀을 가지고 좌표 [math(\textbf{r}_2)]에서 미소부피 [math(d\textbf{r}_2)]안에서 발견될 확률, [math(\cdots)], n번 전자가 [math(m_{sn})] 스핀을 가지고 좌표 [math(\textbf{r}_n)]에서 미소부피 [math(d\textbf{r}_n)]안에서 발견될 확률이| [math(\displaystyle \begin{aligned} |\psi(\textbf{r}_1, \textbf{r}_2, \cdots, \textbf{r}_n, m_{s1}, m_{s2}, \cdots, m_{sn})|^2 d\textbf{r}_1 d\textbf{r}_2 \cdots d\textbf{r}_n \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \sum\limits_{\textup{all} \, m_{s}} |\psi(\textbf{r}_1, \textbf{r}_2, \cdots, \textbf{r}_n, m_{s1}, m_{s2}, \cdots, m_{sn})|^2 d\textbf{r}_1 d\textbf{r}_2 \cdots d\textbf{r}_n \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \Bigg[ \sum\limits_{\textup{all} \, m_{s}}\int \cdots \int \int |\psi(\textbf{r}, \textbf{r}_2, \cdots, \textbf{r}_n, m_{s1}, m_{s2}, \cdots, m_{sn})|^2 d\textbf{r}_2 \cdots d\textbf{r}_n \Bigg] d\textbf{r} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \Bigg[ \sum\limits_{\textup{all} \, m_{s}}\int \cdots \int \int |\psi(\textbf{r}_1, \textbf{r}, \cdots, \textbf{r}_n, m_{s1}, m_{s2}, \cdots, m_{sn})|^2 d\textbf{r}_1 \cdots d\textbf{r}_n \Bigg] d\textbf{r} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \rho(x, y, z) = n\Bigg[ \sum\limits_{\textup{all} \, m_{s}}\int \cdots \int \int |\psi(\textbf{r}, \textbf{r}_2, \cdots, \textbf{r}_n, m_{s1}, m_{s2}, \cdots, m_{sn})|^2 d\textbf{r}_2 \cdots d\textbf{r}_n \Bigg] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \Bigg\langle \psi \Bigg| \sum\limits_{i=1}^n B(\textbf{r}_i) \Bigg| \psi \Bigg\rangle = \sum\limits_{\textup{all} \, m_s} \left[ \int \cdots \int \int \left( \psi^* \sum\limits_{i=1}^n B(\textbf{r}_i) \psi \right) \, d\textbf{r}_1 d\textbf{r}_2 \cdots d\textbf{r}_n \right] = \sum\limits_{i=1}^n \sum\limits_{\textup{all} \, m_s} \left[ \int \cdots \int \int |\psi|^2 B(\textbf{r}_i) \, d\textbf{r}_1 d\textbf{r}_2 \cdots d\textbf{r}_n \right] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \sum\limits_{i=1}^n \sum\limits_{\textup{all} \, m_s} \left[ \int \cdots \int \int |\psi|^2 B(\textbf{r}_i) \, d\textbf{r}_1 d\textbf{r}_2 \cdots d\textbf{r}_n \right] = n \int B(\textbf{r}_1) \sum\limits_{\textup{all} \, m_s} \left[ \int \cdots \int |\psi|^2 \, d\textbf{r}_2 \cdots d\textbf{r}_n \right] \, d\textbf{r}_1 = \int \rho(\textbf{r}_1) B(\textbf{r}_1) \, d \textbf{r}_1 \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \Bigg\langle \psi \Bigg| \sum\limits_{i=1}^n B(\textbf{r}_i) \Bigg| \psi \Bigg\rangle = \int \rho(\textbf{r}) B(\textbf{r}) \, d \textbf{r} \end{aligned} )] |
3. 호헨베르크-콘 이론
3.1. 제1정리
제1정리에선 바닥상태에서 축퇴가 없는 분자의 경우 바닥상태 전자 밀도 [math(\rho_0(\textbf{r}))][4]를 안다면 파동 함수와 모든 물리적 특징들이 하나로 결정됨을 보인다.[5] 이는 "[math(\rho_0(\textbf{r}))]이 해밀토니언을 결정한다."를 보이면 끝나는데 해밀토니언이 결정되면 파동함수와 에너지가 결정되고 파동함수를 통해 물리량 또한 결정되기 때문이다. 먼저 분자의 해밀토니언을 살펴보자. 보른-오펜하이머 근사를 적용시키면 전자의 해밀토니언은| [math(\displaystyle \begin{aligned} \mathcal{H}_e = -\frac{1}{2}\sum\limits_{i=1}^n \nabla_i^2+\sum\limits_{i=1}^n v(\textbf{r}_i)+\sum\limits_{i>j} \frac{1}{r_{ij}} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} v(\textbf{r}_i)=-\sum\limits_{\alpha} \frac{Z_{\alpha}}{r_{i\alpha}} \end{aligned} )] [6] |
| [math(\displaystyle \begin{aligned} E_{0, a} < \langle \psi_{0, b} | \mathcal{H}_a | \psi_{0, b} \rangle &= \langle \psi_{0, b} | \mathcal{H}_a - \mathcal{H}_b + \mathcal{H}_b | \psi_{0, b} \rangle = \langle \psi_{0, b} | \mathcal{H}_a - \mathcal{H}_b | \psi_{0, b} \rangle + E_{0, b} \\ &= \Bigg\langle \psi_{0, b} \Bigg| \sum\limits_{i=1}^n [v_a(\textbf{r}_i)-v_b(\textbf{r}_i ) ] \Bigg| \psi_{0, b} \Bigg\rangle + E_{0, b} \\ &= \int \rho_{0, b}(\textbf{r})[v_a(\textbf{r})-v_b(\textbf{r}) ] d\textbf{r} + E_{0, b} \end{aligned})] |
| [math(\displaystyle \begin{aligned} E_{0, b} < \langle \psi_{0, a} | \mathcal{H}_b | \psi_{0, a} \rangle \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_{0, a} < \int \rho_{0, b}(\textbf{r})[v_a(\textbf{r})-v_b(\textbf{r}) ] d\textbf{r} + E_{0, b} \\ E_{0, b} < \int \rho_{0, a}(\textbf{r})[v_b(\textbf{r})-v_a(\textbf{r}) ] d\textbf{r} + E_{0, a} \end{aligned})] |
| [math(\displaystyle \begin{aligned} E_{0, a} + E_{0, b} < E_{0, a} + E_{0, b} \end{aligned})] |
이제 우리가 가장 큰 관심이 있는 에너지로 가보자. 위 정리에 의해 에너지 또한 밀도의 범함수이다. 이를 [math(E_0 = E_v[\rho_0])][8]로 쓸 것이다.
| [math(\displaystyle \begin{aligned} \mathcal{H}_e = \hat{T} + \hat{V}_{Ne} + \hat{V}_{ee} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_0 = E_v[\rho_0] = \overline{T}[\rho_0]+ \overline{V}_{Ne}[\rho_0]+ \overline{V}_{ee}[\rho_0] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \overline{V}_{Ne} = \Bigg\langle \psi_0 \Bigg| \sum\limits_{i=1}^n v(\textbf{r}_i) ] \Bigg| \psi_0 \Bigg\rangle = \int \rho_0(\textbf{r})v(\textbf{r}) d\textbf{r}\end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_0 = E_v[\rho_0] = \overline{T}[\rho_0]+\overline{V}_{ee}[\rho_0] + \int \rho_0(\textbf{r})v(\textbf{r}) d\textbf{r} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} F[\rho_0] \equiv \overline{T}[\rho_0]+\overline{V}_{ee}[\rho_0] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_0 = E_v[\rho_0] = F[\rho_0] + \int \rho_0(\textbf{r})v(\textbf{r}) d\textbf{r} \end{aligned} )] |
3.2. 제2정리
제2정리는 DFT에서의 변분 원리이다. 다음과 같은 두 조건| [math(\displaystyle \begin{aligned} \int \rho_{\textup{trial}}(\textbf{r}) d\textbf{r} = n \qquad \rho_{\textup{trial}}(\textbf{r}) \geq 0 \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_v[\rho_{\textup{trial}}] \geq E_0 = E_v[\rho_{0}]\end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_v[\rho_{\textup{trial}}] = \langle \psi_{\sf{trial}}|\mathcal{H}_e|\psi_{\sf{trial}} \rangle \geq E_{0} = E_v[\rho_0] \end{aligned} )] |
호헨베르크와 콘이 위 이론들을 증명할 땐 바닥상태의 축퇴도가 없음을 가정했지만 축퇴가 존재하여도 이론을 적용시키는 것이 가능하다.[10]
4. 콘-샴 방법
호헨베르크-콘 이론을 통해 [math(E_0)]를 조건 [math(\int \rho \, d\textbf{r} - N = 0)] 하에서 [math(\rho)]에 대해 최소화 한다면 [math(\rho)]에 대한 방정식을 얻을 것이고 이를 풀면 이론적으론 [math(\rho)]를 파동함수 없이 구할 수 있다. 하지만 [math(\overline{T}[\rho])]와 [math(\overline{V}_{ee}[\rho])]의 꼴을 모르기에[11] 다른 방법을 고안해야 한다. 이 과정에서 핵심 역할을 하는 것이 콘-샴 방법(Kohn-Sham Method)이며 이론상으로 이 방법을 이용하면 정확한 결과값을 얻어낼 수 있다.[12][math(n)]개의 상호작용하지 않는 전자로 이루어져 있으며 모든 전자가 [math(v_s(\textbf{r}_i))]라는 외부 퍼텐셜 하에서 운동하고 있는 가상의 계를 생각해보자(reference system이라고 한다). 여기서 [math(v_s(\textbf{r}_i))]는 reference system의 바닥상태 전자 밀도 [math(\rho_s(\textbf{r}))]가 우리가 관심있어하는 계의 정확한 바닥상태 확률 밀도 [math(\rho_0(\textbf{r}))]와 동일하게 만들어주는 외부 퍼텐셜이다. 호헨베르크-콘 이론에서 증명했듯이 [math(\rho_s(\textbf{r}) = \rho_0(\textbf{r}))]로 전자 밀도를 설정하면 [math(v_s(\textbf{r}_i))]는 하나로 결정된다. 그럼 reference system의 해밀토니언은
| [math(\displaystyle \begin{aligned} \mathcal{H}_s = \sum\limits_{i=1}^n \Big[ -\frac{1}{2}\nabla_i^2 + v_s(\textbf{r}_i) \, \Big] \equiv \sum\limits_{i=1}^n h_i^{KS} \qquad h_i^{KS} \equiv -\frac{1}{2}\nabla_i^2 + v_s(\textbf{r}_i) \end{aligned} )] |
| [math(\displaystyle \begin{aligned} |\psi_{s, 0} \rangle = |u_1^{KS}u_2^{KS} \cdots u_n^{KS} \rangle \end{aligned} )] |
| [math(\displaystyle \begin{aligned} u_i^{KS} = \theta_i^{KS}(\textbf{r}_i)\sigma_i \qquad h_i^{KS} \theta_i^{KS} = \varepsilon_i^{KS} \theta_i^{KS} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} T_s = \sum\limits_{i} \langle \theta^{KS}_i | -\frac{1}{2} \nabla ^ 2 | \theta^{KS}_i \rangle \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_v[\rho_0] = \overline{T}[\rho_0] + \overline{V}_{ee}[\rho_0] + \int \rho_0(\textbf{r})v(\textbf{r}) d\textbf{r} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \Delta T[\rho] \equiv \overline T[\rho] - T_s[\rho] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \Delta V_{ee}[\rho] \equiv \overline{V}_{ee}[\rho] - \frac{1}{2} \int \int \frac{\rho({\textbf{r}}_1)\rho({\textbf{r}}_2)}{r_{12}} d\textbf{r}_1d\textbf{r}_2 = \overline{V}_{ee}[\rho] - J[\rho] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_v[\rho] = T_s[\rho] + J[\rho] + \int \rho(\textbf{r})v(\textbf{r}) d\textbf{r} + \Delta \overline{T}[\rho] + \Delta \overline{V}_{ee}[\rho] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_{xc}[\rho] \equiv \Delta \overline{T}[\rho] + \Delta \overline{V}_{ee}[\rho] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_v[\rho] = T_s[\rho] + J[\rho] + \int \rho(\textbf{r})v(\textbf{r}) d\textbf{r} + E_{xc}[\rho] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \Omega [\{ \theta_i^{KS} \} ] = E_v[\{ \theta_i^{KS} \}] - \sum\limits_{i, j} \varepsilon_{ij} ( \langle i | j \rangle - \delta_{ij} ) \end{aligned})] |
| [math(\displaystyle \begin{aligned} E_v = \sum\limits_{i} \langle i | T | i \rangle + \sum\limits_{i} \langle i | v | i \rangle + \frac{1}{2} \sum\limits_{i, j} \langle ij | ij \rangle + E_{xc} \end{aligned})] |
| [math(\displaystyle \begin{aligned} \Omega = \sum\limits_{i} \langle i | T | i \rangle + \sum\limits_{i} \langle i | v | i \rangle + \frac{1}{2} \sum\limits_{i, j} \langle ij | ij \rangle + E_{xc} - \sum\limits_{i, j} \varepsilon_{ij} ( \langle i | j \rangle - \delta_{ij} ) \end{aligned})] |
| [math(\displaystyle \begin{aligned} \delta \Omega &= \sum\limits_{i} ( \langle \delta i | T | i \rangle + \langle i | T | \delta i \rangle) + \sum\limits_{i} ( \langle \delta i | v | i \rangle + \langle i | v | \delta i \rangle ) \\ &+ \sum\limits_{i, j} ( \langle \delta ij | ij \rangle + \langle ij | \delta ij \rangle ) + \delta E_{xc} - \sum\limits_{i, j} ( \varepsilon_{ij} \langle \delta i | j \rangle - \varepsilon_{ji} \langle j | \delta i \rangle ) = 0 \end{aligned})] |
| [math(\displaystyle \begin{aligned} \delta E_{xc} = \int \frac{\delta E_{xc}}{\delta \rho} \delta \rho \, d \textbf{r} \end{aligned})] |
| [math(\displaystyle \begin{aligned} \rho(\textbf{r}) = \sum\limits_{i} \theta^{KS*}_i(\textbf{r})\theta^{KS}_i(\textbf{r}) \end{aligned})] |
| [math(\displaystyle \begin{aligned} \delta \rho = \sum\limits_{i} \theta^{KS*}_i \delta \theta^{KS}_i + \sum\limits_{i} \theta^{KS}_i \delta \theta^{KS*}_i \end{aligned})] |
| [math(\displaystyle \begin{aligned} \delta E_{xc} = \int \frac{\delta E_{xc}}{\delta \rho} \left( \sum\limits_{i} \theta^{KS*}_i \delta \theta^{KS}_i + \sum\limits_{i} \theta^{KS}_i \delta \theta^{KS*}_i \right) \, d \textbf{r} \end{aligned})] |
| [math(\displaystyle \begin{aligned} v_{xc} \equiv \frac{\delta E_{xc}}{\delta \rho} \end{aligned})] |
| [math(\displaystyle \begin{aligned} \delta E_{xc} = \sum\limits_{i} ( \langle \delta i | v_{xc} | i \rangle + \langle i | v_{xc} | \delta i \rangle ) \end{aligned})] |
| [math(\displaystyle \begin{aligned} \bigg[-\frac{1}{2}\nabla^2 + v(\textbf{r}) + \int \frac{\rho(\textbf{r}')}{|\textbf{r} - \textbf{r}'|} d \textbf{r}' + v_{xc}(\textbf{r}) \bigg]\theta_i^{KS}(\textbf{r}) = \varepsilon_i^{KS} \theta_i^{KS}(\textbf{r}) \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \left( -\frac{1}{2}\nabla^2 + v_s(\textbf{r}) \right) \theta_i^{KS}(\textbf{r}) = \varepsilon_i^{KS} \theta_i^{KS}(\textbf{r}) \end{aligned} )] |
| [math(\displaystyle \begin{aligned} v_s(\textbf{r}) = v(\textbf{r}) + \int \frac{\rho(\textbf{r}')}{|\textbf{r} - \textbf{r}'|} d \textbf{r}' + v_{xc}(\textbf{r}) = - \sum\limits_{\alpha} \frac{Z_{\alpha}}{r_{i\alpha}} + \int \frac{\rho(\textbf{r}')}{|\textbf{r} - \textbf{r}'|} d \textbf{r}' + v_{xc}(\textbf{r}) \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \left( -\frac{1}{2}\nabla^2 + v_\textup{eff}(\textbf{r}) \right) \theta_i^{KS}(\textbf{r}) = \varepsilon_i^{KS} \theta_i^{KS}(\textbf{r}) \end{aligned} )] |
- [math(\rho(\textbf{r}))]을 가정한다. 보통 원자의 전자 밀도를 중첩시켜서 얻는다. 가정한 [math(\rho(\textbf{r}))]를 통해 [math(v(\textbf{r}))]와 [math(v_{xc}(\textbf{r}))]를 계산한다.
- 위에서 얻은 전자밀도를 방정식에 넣은 뒤 푼다. 방정식을 수치적으로 풀거나[16] 하트리-포크 방법에서 사용한 basis set approach를 사용하여 basis들의 계수를 찾는다.
- [math(\rho(\textbf{r}) = \sum\limits_{i} |\theta^{KS}_i (\textbf{r})|^2)]을 이용하여 [math(\rho(\textbf{r}))]을 재추정한다.
- [math(\rho(\textbf{r}))]와 [math(E_{xc}[\rho])]가 수렴할 때 까지 1~3번 과정을 반복한다.
5. 교환-상관 에너지
콘-샴 방정식의 가장 큰 문제점은 [math(E_{xc})]와 [math(v_{xc})]의 모양을 모른다는 것이다. DFT 정확도의 핵심은 이 퍼텐셜을 실험치와 잘 부합하도록 근사하는 것이다. 그렇기에 정말 수많은 종류의 [math(E_{xc})]와 [math(v_{xc})]가 존재하며[17] 지금까지도 연구되고 있다. 아래는 그중 자주 쓰이는 일부만을 추려낸 것이다.5.1. LDA
가장 기초적인 방법이 균일한 전자 가스로 부터 [math(E_{xc})]를 근사하는 것이다. 이를 LDA(Local Density Approximation)라 한다. 여기서 [math(E_{xc})]는| [math(\displaystyle \begin{aligned} E_{xc}^{LDA}[\rho] = \int \rho(\textbf{r}) \varepsilon_{xc}(\rho) d\textbf{r} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} v_{xc}^{LDA} = \frac{\delta E_{xc}^{LDA}}{\delta \rho} = \varepsilon_{xc}(\rho) + \rho \frac{\partial \varepsilon_{xc}}{\partial \rho(\rho)} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \varepsilon_{xc}(\rho) = \varepsilon_{x}(\rho) + \varepsilon_{c}(\rho) \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_{xc} = E_{x} + E_{c} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_{x}[\rho] &= \int \rho(\textbf{r}) \varepsilon_{x}(\rho) d\textbf{r} \\ E_{c}[\rho] &= \int \rho(\textbf{r}) \varepsilon_{c}(\rho) d\textbf{r}\end{aligned} )] |
5.2. LSDA
open-shell molecule들이나 해리 상태에 가까운 molecule들의 경우 LSDA(Local Spin Density Approximation)을 사용하는 것이 LDA보다 더 좋은 결과를 준다. HF의 UHF처럼 같은 콘-샴 오비탈에 들어가더라도 다른 모양을 줌으로써 오비탈의 유연성을 늘린다. [math(\alpha)] 스핀을 가지는 전자들의 밀도를 [math(\rho^{\alpha}(\textbf{r}))], [math(\beta)] 스핀을 가지는 전자들의 밀도를 [math(\rho^{\beta}(\textbf{r}))]라 하면 [math(E_{xc})]는 [math(\rho^{\alpha}(\textbf{r}))]와 [math(\rho^{\beta}(\textbf{r}))] 두개의 함수의 범함수가 된다.| [math(\displaystyle \begin{aligned} E_{xc} = E_{xc}[\rho^{\alpha}, \rho^{\beta}] \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \rho = \rho^{\alpha} + \rho^{\beta} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} \zeta = \frac{\rho^{\alpha}-\rho^{\beta}}{\rho^{\alpha}+\rho^{\beta}} \end{aligned} )] |
- Slater-Dirac exchange functional(Slater)[19]
[math(\displaystyle \begin{aligned} E^{\textup{Slater}}_x[\rho^{\alpha}, \rho^{\beta}] &= \int \rho \varepsilon_x(\rho, \zeta) d \textbf{r} \\ \varepsilon_x(\rho, \zeta) &= \varepsilon^0_x(\rho) + [\varepsilon^1_x(\rho) -\varepsilon^0_x(\rho) ] f(\zeta) \\ \varepsilon^0_x(\rho) &= \varepsilon_x(\rho, 0) = C_x \rho^{1/3} \\ \varepsilon^1_x(\rho) &= \varepsilon_x(\rho, 1) = 2^{1/3}C_x \rho^{1/3} \\ C_x &= \frac{3}{4}\left( \frac{3}{\pi} \right)^{1/3} \\ f(\zeta) &= \frac{(1+\zeta)^{4/3}+(1-\zeta)^{4/3}-2}{2(2^{1/3}-1)} \\ v^{\textup{Slater}}_{x\sigma} &= \frac{\delta E^{\textup{Slater}}_{x}}{\delta \rho_{\sigma}} = \left( \frac{6}{\pi} \rho_{\sigma} \right) ^ {1/3} \end{aligned} )] - Vosko-Wilk-Nusair correlation functional 5# (VWN5)[20]
[math(\displaystyle \begin{aligned} E^{\textup{VWN}}_c[\rho^{\alpha}, \rho^{\beta}] &= \int \rho \varepsilon_c^{\textup{VWN}}(\rho, \zeta) d \textbf{r} \\ \varepsilon_c^{\textup{VWN}}(\rho, \zeta) &= \varepsilon_\textup{I}(r_s) + \Delta \varepsilon_c(r_s, \zeta) \\ \varepsilon_i(r_s) &= A_i \left[ \textup{ln} \frac{x^2}{X(x)} +\frac{2b}{Q} \textup{tan}^{-1}\left( \frac{Q}{2x+b} \right) - \frac{bx_0}{X(x_0)} \left\{ \textup{ln} \frac{(x-x_0)^2}{X(x)} + \frac{2(b+2x_0)}{Q} \textup{tan}^{-1} \left( \frac{Q}{2x+b} \right) \right\} \right] \\ x &= r_s ^ {1/2} \\ \frac{4}{3}\pi r_s ^3 &= \frac{1}{\rho} \\ Q &= (4c_i - b_i^2)^{1/2} \\ X(x) &= x^2 + b_ix + c_i(i= \textup{I}, \textup{II} ) \\ \Delta \varepsilon_c(r_s, \zeta) &= \varepsilon_{\textup{III}}(r_s) \left[ \frac{f(\zeta)}{f(0)} \right] \left[ 1+\beta_i(r_s) \zeta ^ 4 \right] \\ \beta_i(r_s) &= \left[ \frac{f(0)}{\varepsilon_{\textup{III}}(r_s)} \right] \Delta \varepsilon(r_s, 1) - 1 \\ \Delta \varepsilon(r_s, 1) &= \varepsilon_\textup{I}(r_s) - \varepsilon_\textup{II}(r_s) \end{aligned} )]
5.3. GGA
LDA와 LSDA는 균일한 전자 가스로 부터 유도 되었기에 [math(\rho)]가 천천히 변할 때만 적절하다. 그렇기에 GGA(Generalized-Gradient Approximation)(Gradient-corrected functional이라고도 불린다.)에선 [math(\nabla \rho)] 항을 추가하여[21] 정확도를 높인다. 이를 식으로 표현하면| [math(\displaystyle \begin{aligned} E^{\textup{GGA}}_{xc}[\rho^{\alpha}, \rho^{\beta}] = \int f(\rho^{\alpha}, \rho^{\beta}, \nabla \rho^{\alpha}, \nabla \rho^{\beta}) d \textbf{r} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E^{\textup{GGA}}_{xc} = E^{\textup{GGA}}_{x} + E^{\textup{GGA}}_{c} \end{aligned} )] |
- Becke exchange functional from 1988 (B88)[22]
[math(\displaystyle \begin{aligned} E^{\textup{B88}}_{x} &= E^{\textup{LSDA}}_x - b \sum\limits_{\sigma = \alpha, \beta} \int \frac{(\rho^{\sigma})^{4/3}\chi^2_{\sigma}}{1+6b\chi_{\sigma}\textup{ln}[\chi_{\sigma} + (\chi^2_{\sigma}+1)^{1/2} ]} d \textbf{r} = E^{\textup{LSDA}}_x + \Delta E^{\textup{B88}}_x \\ \chi_{\sigma} &= \frac{|\nabla \rho^{\sigma}|}{(\rho^{\sigma})^{4/3}} \\ E^{\textup{LSDA}}_x &= -\frac{3}{4} \left( \frac{6}{\pi} \right) ^ {1/3} \int \left[ \left( \rho^{\alpha} \right) ^ {4/3} + \left( \rho^{\beta} \right) ^ {4/3} \right] d \textbf{r} \end{aligned} )]
b는 조정 가능한 매개변수로 원자 단위계에선 0.0042이다. - Perdew-Burke-Ernzerhof exchange functional (PBE)[23], Revised PBE (revPBE)[24]
3\pi^2 \rho \right)^{1/3}\end{aligned} )]||
여기서 쓰이는 매개변수 값은 두 방법 모두[math(\mu = 0.21951)]이고 PBE에선 [math(\kappa = 0.804)], revPBE에선 [math(\kappa = 1.245)]를 사용한다. - Perdew-Burke-Ernzerhof correlation functional (PBE)[25]
[math(\displaystyle \begin{aligned} E^{\textup{PBE}}_{c} &= E^{\textup{LSDA}}_c + \int H^{\textup{PBE}}(r_s, \zeta, t) d \textbf{r} = E^{\textup{LSDA}}_c + \Delta E^{\textup{PBE}}_c \\ H^{\textup{PBE}} &= \gamma \phi^3 \textup{ln} \left\{ 1+ \frac{\beta}{\gamma}t^2 \left( \frac{1+At^2}{1+At^2+ A^2t^4} \right) \right\} \\ A &= \frac{\beta}{\gamma}\{ \textup{exp}(-\varepsilon^{\textup{LSDA}}_c/\gamma \phi ^ 3) -1 \}^{-1} \\ t &= \frac{|\nabla \rho|}{2k_s\phi\rho} \\ \phi &= \frac{\left( 1 + \zeta \right)^{2/3} + \left( 1 - \zeta \right) ^ {2/3}}{2} \\ k_s &= \left( \frac{4k_F}{\pi}\right)^{1/2} \end{aligned} )]
여기서 쓰이는 매개변수 값은 [math(\gamma = 0.031091, \beta = 0.066725)]이다. - Lee-Yang-Parr opposite-spin correlation functional (LYP)[26]
[math(\displaystyle \begin{aligned} E^{\textup{LYP}}_c [\rho^{\alpha}, \rho^{\beta}] &= -a \int \frac{\gamma}{1+d\rho^{-1/3}} \left[ \rho + 2b\rho^{-5/3} \left\{ 2^{2/3} C_F \rho^{8/3}_{\alpha} + 2^{2/3} C_F \rho^{8/3}_{\beta} -pt_w + \frac{1}{9} \left( \rho_{\alpha}t^{\alpha}_w + \rho_{\beta}t^{\beta}_w \right) + \frac{1}{18} \left( \rho_{\alpha} \nabla^2 \rho_{\alpha} + \rho_{\beta} \nabla^2 \rho_{\beta} \right) \right\} \textup{exp} \left( -c\rho^{-1/3} \right) \right] \\ \gamma &= 2 \left( 1-\frac{\rho^2_{\alpha}+ \rho^2_{\beta}}{\rho^2} \right) \\ t_w &= \frac{1}{8}\frac{|\nabla \rho|^2}{\rho} - \frac{1}{8}\nabla^2 \rho \\ C_F &= \frac{3}{10} \left( 3\pi^2 \right)^{2/3} \end{aligned} )]
[math(a = 0.04918, b = 0.132, c = 0.2533, d = 0.349)]
5.4. Meta-GGA
GGA에서 더 정확도를 높이기 위해 meta-GGA에선 [math(\nabla^2 \rho)]와 아래와 같은 운동 에너지 밀도(Kinetic-energy density)를 추가한다.| [math(\displaystyle \begin{aligned} \tau_{\alpha} &= \frac{1}{2} \sum\limits_{i} |\nabla \theta^{\textup{KS}}_{i{\alpha}}|^2 \\ \tau_{\beta} &= \frac{1}{2} \sum\limits_{i} |\nabla \theta^{\textup{KS}}_{i{\beta}}|^2 \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E^{\textup{MGGA}}_{xc}[\rho^{\alpha}, \rho^{\beta}] = \int f(\rho^{\alpha}, \rho^{\beta}, \nabla \rho^{\alpha}, \nabla \rho^{\beta}, \nabla^2 \rho^{\alpha}, \nabla^2 \rho^{\beta}, \tau_{\alpha}, \tau_{\beta}) d \textbf{r} \end{aligned} )] |
- Tao-Perdew-Staroverov-Scuseria exchange, correlation functional (TPSS)[27], Revised TPSS (revTPSS)[28]
- M06-L[29][30], revM06-L
5.5. Hybrid Functional
Hybrid functional에선 [math(E_{xc})]에 하트리-포크 방법의 [math(E_x)]를 추가한다. [math(E_x^{\textup{GGA}})]에 추가하냐 [math(E_x^{\textup{meta-GGA}})]에 추가햐냐에 따라 hybrid GGA와 hybrid meta-GGA로 세분화하기도 한다.20% HF exchange + 8% Slater LSDA exchange + 72% B88 GGA exchange + 19% VWN LSDA correlation + 81% LYP GGA correlation
| [math(\displaystyle \begin{aligned} E^{\textup{B3LYP}}_{xc} &= 0.20E^{\textup{HF}}_{x} + 0.08E^{\textup{Slater}}_{x} + 0.72E^{\textup{B88}}_{x} + 0.19E^{\textup{VWN}}_{c} + 0.81E^{\textup{LYP}}_{c} \\ &= 0.20E^{\textup{HF}}_{x} + 0.80E^{\textup{Slater}}_{x} + 0.72\Delta E^{\textup{B88}}_{x} + 0.19E^{\textup{VWN}}_{c} + 0.81E^{\textup{LYP}}_{c} \end{aligned} )] |
- PBE0[33], revPBE0
25% HF exchange + 75% PBE GGA exchange + PBE GGA correlation
25% HF exchange + 75% revPBE GGA exchange + PBE GGA correlation[math(\displaystyle \begin{aligned} E^{\textup{PBE0}}_{xc} &= 0.25E^{\textup{HF}}_{x} + 0.75E^{\textup{PBE}}_{x} + E^{\textup{PBE}}_{c} \\ E^{\textup{revPBE0}}_{xc} &= 0.25E^{\textup{HF}}_{x} + 0.75E^{\textup{revPBE}}_{x} + E^{\textup{PBE}}_{c} \end{aligned} )]
- B97[34]
- TPSSh[39], revTPSSh[40]
10% HF exchange + 90% TPSS meta-GGA exchange + TPSS meta-GGA correlation
10% HF exchange + 90% revTPSS meta-GGA exchange + revTPSS meta-GGA correlation[math(\displaystyle \begin{aligned} E^{\textup{TPSSh}}_{xc} &= 0.10E^{\textup{HF}}_{x} + 0.90E^{\textup{TPSS}}_{x} + E^{\textup{TPSS}}_{c} \\ E^{\textup{revTPSSh}}_{xc} &= 0.10E^{\textup{HF}}_{x} + 0.90E^{\textup{revTPSS}}_{x} + E^{\textup{revTPSS}}_{c} \end{aligned} )]
5.6. Double-Hybrid Functional
Double Hybrod Functional에선 unoccupied KS orbitals 의존성을 추가해[41] Hybrid Functional보다 정확도를 높인다.- Becke, 2-parameter, Perturbation, Lee–Yang–Parr (B2PLYP)[42]
53% HF exchange + 47% B88 GGA exchange + 73% LYP GGA correlation + 27% MP2 correlation[math(\displaystyle \begin{aligned} E^{\textup{B2PLYP}}_{xc} = 0.53E^{\textup{HF}}_{x} + 0.47E^{\textup{B88}}_{x} + 0.73E^{\textup{LYP}}_{c} + 0.27E^{\textup{MP2}}_{c}\end{aligned} )] - XYG3[43][44]
80.33% HF exchange - 1.4% Slater LSDA exchange + 21.07% B88 GGA exchange + 67.89% LYP GGA correlation + 32.11% MP2 correlation[math(\displaystyle \begin{aligned} E^{\textup{XYG3}}_{xc} = 0.8033E^{\textup{HF}}_{x} - 0.014E^{\textup{Slater}}_{x} + 0.2107E^{\textup{B88}}_{x} + 0.6789E^{\textup{LYP}}_{c} + 0.3211E^{\textup{MP2}}_{c}\end{aligned} )] - XYGJ-OS[45]
- DSD-PBEPBE-D3[46]
- PTPSS-D3[47]
6. 한계
DFT에도 문제점들이 있는데 그중 하나가 장거리 전자 상관(Long-Range Electronic Correlation)을 고려하지 못한다는 것이다. 단거리 전자 상관(Short-Range Electronic Correlation)을 표현하는 것에만 신경을 쓰므로 공유 결합과 같은 것들은 잘 설명할지 몰라도 장거리 전자 상관이 중요한 전하 이동(Charge Transfer)나 분산력(Dispersion Force) 등을 설명하지 못한다는 단점이 있다. 이를 해결하기 위한 방법 중 하나론 RSH(Range-Separated Hybrid)이라 불리는 방법으로 다음과 같이 단거리에선 DFT Exchange를, 장거리에선 HF Exchange를 사용하는 것이다. 대표적인 예시인 CAM-B3LYP[48]에선| [math(\displaystyle \begin{aligned} \frac{1}{r_{12}} = \underbrace{\frac{1 - [\alpha + \beta \textup{erf}(\omega r_{12}) ]}{r_{12}}}_{\textup{Short Range: DFT Exchange}} + \underbrace{\frac{\alpha + \beta \textup{erf}(\omega r_{12}) }{r_{12}}}_{\textup{Long Range: HF Exchange}} \end{aligned} )] |
| [math(\displaystyle \begin{aligned} E_{x} = E_{x}^{SR} + E_{x}^{LR} = \frac{1 - [\alpha + \beta \textup{erf}(\omega r_{12}) ]}{r_{12}} E_x^{\textup{B88}} + \frac{\alpha + \beta \textup{erf}(\omega r_{12}) }{r_{12}} E_x^{\textup{HF}}\end{aligned} )] |
- wB97X-V(RSH GGA)[49]
- wB97X-D3(RSH GGA)[50]
- wB97X-D(RSH GGA)[51]
- wB97M-V(RSH meta-GGA)[52]
또 다른 방법인 DFT-D(Dispersion-corrected DFT)에선 DFT로 구한 에너지에 실험적으로 구한 분산력 보정치를 더해준다. 예시로 DFT-D2에선 아래와 같은 보정항을 추가한다.
| [math(\displaystyle \begin{aligned} E_{\textup{DFT-D2}} = E_{\textup{DFT}} - s_6 \sum\limits_{A<B} \frac{C^{AB}_6}{|\textbf{R}_A - \textbf{R}_B|^6} f^{\textup{D2}}_{dmp}(|\textbf{R}_A - \textbf{R}_B|) \\ f^{\textup{D2}}_{dmp}(|\textbf{R}_A - \textbf{R}_B|) = \left[ 1 + e ^ {-d \left( {|\textbf{R}_A - \textbf{R}_B| / R_{0, AB} - 1} \right) } \right]^{-1} \end{aligned} )] |
7. 하트리-포크 방법과의 비교
- 근사 유무
하트리-포크 방정식은 오비탈 근사를 도입했기 때문에 방정식을 정확하게 풀더라도 항상 분자의 실제 에너지보다는 높게 나온다(HF limit) 하지만 콘-샴 방정식의 경우 [math(v_{xc})]만 정확하다면 정확한 에너지를 내어준다.[53] - 개선의 방향성
하트리-포크 방법의 경우 들뜬 상태의 Slater Determinant를 추가해 줌으로써 정확한 해에 다가갈 수 있다. 하지만 DFT의 경우 [math(v_{xc})]를 개선하는 정해진 가이드라인이 없다. - 변분 원리
하트리-포크 방정식의 경우 오비탈 근사 이외에는 사용한 근사가 없기에 변분 원리를 따른다. DFT 또한 제2정리에 의해 변분 원리를 따르지만 [math(v_{xc})]를 근사해서 넣기에 구한 에너지가 참값보다 작게 나올수도 있다.
8. 관련 문서
[1] 공간에 관한 좌표[math(3N)]개, 스핀에 관한 좌표 [math(N)]개[2] 핵의 위치, 핵의 전하량[3] 이는 [math(|\psi|^2)]가 두 입자의 변환에 대해 불변임을 이용하면 쉽게 증명 가능하다.[4] 0은 바닥상태임을 나타낸다.[5] 즉 앞서 언급한 것들은 [math(\rho(\textbf{r}))]의 범함수이다.[6] 이 양을 DFT에선 외부 퍼텐셜(external potential)이라 한다.[7] 상수 차이만 나는 퍼텐셜은 같은 퍼텐셜로 본다.[8] [math(v)]는 [math(E_0)]가 [math(v(\textbf{r}))]에 의존한다는 것을 강조하기 위해 붙었다.[9] 사실 이것은 [math(\rho_{\textup{trial}})]에 해당하는 반대칭 파동함수를 주는 [math(v_{\textup{trial}})]가 존재할때만 성립한다. 이러한 [math(\rho_{\textup{trial}})]를 [math(v)]-representable하다고 하며 모든 [math(\rho_{\textup{trial}})]가 [math(v)]-representable하지는 않다. 다만 크게 걱정할 필요는 없는게 [math(v)]-representable 조건을 약화시킬 수 있는데 이는 Density-Finctional Theory of Atoms and Molecules를 참고하여라.[10] 이 또한 Density-Finctional Theory of Atoms and Molecules를 참고하여라.[11] 물론 이 항을 근사적으로 나타낸다면 이 방법 또한 [math(\rho)]를 구하는데 쓸 수 있다. 이를 Orbital-free density functional theory라 한다. 하지만 정확도 문제로 자주 쓰이는 방법은 아니다.[12] 물론 이론상으로 그렇다는 것이고 실제로는 근사가 필요하다.[13] 이제부터 첨자 0은 빼도록 하겠다.[14] 전기 퍼텐셜 참조[15] 아래에서 나오는 표기법들은 하트리-포크 방법 참고[16] 잘 쓰이는 방법은 아니다.[17] 참고[18] 오일러-라그랑주 방정식 참조[19] P. A. M. Dirac (1930) P. Camb. Philos. Soc. 26, pp. 376.[20] S. H. Vosko, L. Wilk, and M. Nusair (1980) Can. J. Phys. 58, pp. 1200.[21] [math(\nabla^2 \rho)]까지 추가하기도 한다.[22] A. D. Becke (1988) Phys. Rev. A 38, pp. 3098.[23] J. P. Perdew, K. Burke, and M. Ernzerhof (1996) Phys. Rev. Lett. 77, pp. 3865.[24] Y. Zhang and W. Yang (1998) Phys. Rev. Lett. 80, pp. 890.[25] J. P. Perdew, K. Burke, and M. Ernzerhof (1996) Phys. Rev. Lett. 77, pp. 3865.[26] C. Lee, W. Yang, and R. G. Parr (1988) Phys. Rev. B 37, pp. 785.[27] J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria (2003) Phys. Rev. Lett. 91, pp. 146401.[28] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Constantin, and J. Sun (2009) Phys. Rev. Lett. 103, pp. 026403.[29] Y. Zhao and D. G. Truhlar (2006) J. Chem. Phys. 125, pp. 194101.[30] 전이 금속, 무기화합물, 유기금속화합물 계산에 주로 쓰인다.[31] A. D. Becke (1993) J. Chem. Phys. 98, pp. 5648.[32] P. J. Stephens, F. J. Devlin, C. F. Chabolowski, and M. J. Frisch (1994) J. Phys. Chem. 98, pp. 11623.[33] C. Adamo and V. Barone (1999) J. Chem. Phys. 110, pp. 6158.[34] A. D. Becke (1997) J. Chem. Phys. 107, pp. 8554.[35] Y. Zhao and D. G. Truhlar (2008) Theor. Chem. Acc. 120, pp. 215.[36] 주족 원소 열화학, 동역학, 비공유결합성 상호작용 계산에선 매우 좋은 성능을 지니지만 전이 금속 열화학이나 유기금속화합물 계산엔 사용할 수 없다.[37] Y. Zhao and D. G. Truhlar (2007) J. Chem. Theory Comput. 4, pp. 1849.[38] 주족 원소 열화학, 동역학, 비공유결합성 상호작용 계산에 쓰인다.[39] V. N. Staroverov, G. E. Scuseria, J. Tao, and J. P. Perdew (2003) J. Chem. Phys. 119, pp. 12129.[40] G. I. Csonka, J. P. Perdew, and A. Ruzsinszky (2010) J. Chem. Theory Comput. 6, pp. 3688.[41] 한 예시로 MP2의 이차 보정에너지를 추가한다.[42] S. Grimme (2006) J. Chem. Phys. 124, pp. 034108.[43] Y. Zhang, X. Xu, and W. A. Goddard III (2009) Proc. Natl. Acad. Sci. USA 106, pp. 4963.[44] 특이하게 [math(E^{\textup{MP2}}_{c})]의 계산엔 B3LYP의 오비탈을 사용한다.[45] I. Y. Zhang, X. Xin, Y. Jung, and W. A. Goddard III (2011) Proc. Natl. Acad. Sci. USA 108, pp. 19896.[46] S. Kozuch and J. M. L. Martin (2013) J. Comput. Chem. 34, pp. 2327.[47] L. Goerigk and S. Grimme (2011) J. Chem. Theory Comput. 7, pp. 291.[48] T. Yanai, D. P. Tew, and N. C. Handy (2004) Chem. Phys. Lett. 393, pp. 51.[49] N. Mardirossian and M. Head-Gordon (2014) Phys. Chem. Chem. Phys. 16, pp. 9904.[50] Y.-S. Lin, G.-D. Li, S.-P. Mao, and J.-D. Chai (2013) J. Chem. Theory Comput. 9, pp. 263.[51] J.-D. Chai and M. Head-Gordon (2008b) Phys. Chem. Chem. Phys. 10, pp. 6615.[52] N. Mardirossian and M. Head-Gordon (2016) J. Chem. Phys. 144, pp. 214110.[53] 이론상 그렇다는 것이고 실제론 [math(v_{xc})]를 정확하게 찾는 것은 불가능하다.