연구 관련 자료

제일 원리적 전자 구조 계산법

고체계의 여러 물리적 성질은 그 계의 총에너지와 어떤 파라미터에 대한 그 미분값과 밀접한 관련이 있다. 예를 들어 절대온도 0K에서 고체나 표면은 총에너지가 가장 낮은 구조를 가지며, 평형 구조 근처에서 격자 상수의 길이에 대한 미분으로부터 bulk modulus, 포논 주기등을 계산할 수 있다. 그래서, 어떤 계의 총에너지를 정확히 계산하는 것은 중요한 문제가 된다.

총에너지를 계산하기 위해서는 고체를 구성하는 이온과 전자들간의 상호작용을 정확히 기술할 수 있어야 한다. 우리가 사용하는 방법의 이론적 근거는 밀도 범함수 이론 (Density Functional Theory), 슈도퍼텐셜, 그리고 역격자 공간에서의 총에너지 전개방법으로써 항목별로 아래에 상술하겠다.

1. 밀도 범함수 이론 (Density Functional Theory)

총에너지계산에서 밀도 범함수 이론은 큰 역할을 한다. 이것은 상호 작용하는 전자계의 기저 상태 에너지가 전자 밀도에 의해 유일하게 결정된다는 것으로, 그 시발점을 1920년대 Thomas와 Fermi로 부터 찾을 수 있다. 그러나, 본격적으로 연구가 된 것은 1964년 Hohenberg 와 Kohn 의 논문[1]이 발표된 후이다. 그 후 DFT는 원자나 분자 그리고 고체의 성질을 연구할 때 중요한 이론 근거가 되어 왔으며, 1965년에는 Kohn과 Sham[2]이 전자의 다체문제를 자체 충족적인 단일 전자 문제로 바꿀 수 있음을 보여주어 실제적인 고체 계산의 길이 열었다.

N개의 전자들이 상호작용하는 계의 해밀토니안은 다음과 같이 주어진다.
 

(1)


 

여기에서 $ T = \sum_{i=1}^N (-\frac12 \nabla_i^2) $ 로 운동에너지를, $V_{\rm ext}$ 는 외부 퍼텐셜을, 그리고 $V_{\rm ee} = \sum_{i<j}^N \frac{1}{r_{ij}}$ 으로 전자상호작용에너지를 나타낸다. 우리는 여기서 $\hbar = m = e = 1$인 원자 단위를 사용한다.

우리가 위와 같은 계를 이해하기 위해서는 다체계의 슈뢰딩거 방정식을 풀어 복잡한 N-전자파동함수 $\Psi({\bf r}_1,{\bf r}_2,\cdots,{\bf r}_N)$ 을 구해야 하지만, DFT는 파동함수를 직접 구하는 것이 아니라, 훨씬 더 간단한 양인 전자밀도로 계의 성질을 나타낼 수 있다는 것이다.

Hohenberg와 Kohn은 다음과 같은 두가지 정리를 발표하였다. 그 첫번째 정리는 전자밀도를 이용하여 전자가 상호작용하는 계의 기저 상태를 기술할 수 있다는 것이다. 즉 어떤 해밀토니안이 주어졌을때, 외부퍼텐셜은 전자밀도에 의해서 상수를 제외하고 결정이 된다는 것이다. 그리고, 그 밀도에 해당되는 총에너지 $E_v[\rho]$는 다음과 같이 쓸 수 있다.
 

\begin{displaymath} (2)


 

여기서
 

\begin{displaymath}
F[\rho] = \frac12\int d{\bf r}\int d{\bf r'}
\frac{\rho({\...
...)}{\vert{\bf r}- {\bf r'}\vert}
+ T[\rho] + E_{\rm xc}[\rho]
\end{displaymath} (3)


 

식 (3)에서, 첫째항은 고전적인 쿨롱 에너지(하트리 에너지)를, $T[\rho]$, $E_{\rm xc}[\rho]$ 는 각각 운동에너지, 교환-상관에너지를 나타낸다. 한편, $F[\rho]$$V_{\rm ext}({\bf r})$에 관계가 없는 보편적인 함수가 된다.

Hohenberg-Kohn의 두번째 정리는 에너지 변분원리로 다음과 같다. 즉, $\rho'({\bf r}) \ge 0$$\int \rho'({\bf r}) d {\bf r} = N$을 만족하는 전자밀도 $\rho'({\bf r})$ 에 대하여 다음이 성립한다.
 

\begin{displaymath}
E_0 \le E_v[\rho']
\end{displaymath} (4)


 

여기서 $E_0$는 기저 상태 에너지이다. 그리고 식 (2)에 변분원리를 적용하면 다음과 같은 Euler 식을 얻는다.
 

\begin{displaymath}
\mu = V_{\rm ext}({\bf r}) + \frac{\delta F[\rho]}{\delta \rho({\bf r})},
\end{displaymath} (5)


 

여기서 $\mu$는 전자의 수가 일정하다는 구속조건을 만족시키기 위한 라그랑지안 계수이다. 우리는 식(5)를 풀어서 에너지를 최소화하는 밀도를 구해야한다.

한편, Kohn과 Sham은 운동에너지항을 상호작용하지 않는 계의 운동에너지항 $T_{\rm s}[\rho]$로 근사하고, 둘의 차이를 $E_{\rm xc}$에 포함시켜, 다음과 같은 단일 전자 방정식을 유도하였다. 이것을 Kohn-Sham 방정식이라 부른다.
 

\begin{displaymath}
\left [-\frac12 \nabla^2 + V_{\rm eff}({\bf r})\right] \psi_i({\bf r})
= \varepsilon_i \psi_i({\bf r}),
\end{displaymath} (6)


 

그리고 밀도는 다음과 같이 주어진다.
 

\begin{displaymath}
\rho({\bf r})= \sum_i^N \vert\psi_i({\bf r})\vert^2.
\end{displaymath} (7)


 

식 (6)에서 $V_{\rm eff}({\bf r})= V_{\rm ext}({\bf r})+ V_{\rm H}({\bf r})+V_{\rm xc}({\bf r})$ 로 유효퍼텐셜을 나타낸다. 그리고, 하트리 퍼텐셜($V_{\rm H}$)와 교환-상관퍼텐셜($V_{\rm xc}$)는 다음과 같다.
 

$\displaystyle V_{\rm H}({\bf r})$ $\textstyle =$ $\displaystyle \int d {\bf r'} \frac{\rho({\bf r'})}{\vert{\bf r}- {\bf r'}\vert},$ (8)
$\displaystyle V_{\rm xc}({\bf r})$ $\textstyle =$ $\displaystyle \frac{\delta E_{\rm xc}[\rho]}{\delta \rho({\bf r})}.$ (9)

위의 식 (6)과 (7)에서 $V_{\rm eff}$$\rho({\bf r})$는 서로 의존하기 때문에 자체 충족적으로 풀어야 한다. 그래서, 우선 $\rho({\bf r})$을 추정하여 식(6)를 풀고, 새로운 $\rho({\bf r})$을 계산한다. 이 과정을 반복하여, 입력과 출력밀도의 차이가 일정한 기준안에 들어오면 계산을 멈춘다.

한편, 총에너지는 다음과 같이 쓰여진다.
 

$\displaystyle E[\rho]$ $\textstyle =$ $\displaystyle \int d{\bf r}\rho({\bf r})V_{\rm ext}({\bf r})
+ \sum_i^N \int d {\bf r}\psi_i^{\ast}\left( -\frac12 \nabla^2\right)\psi_i$ (10)
  $\textstyle +$ $\displaystyle \frac12 \int d{\bf r}\int d{\bf r'}\frac{\rho({\bf r})\rho({\bf r})}{\vert{\bf r}-{\bf r'}\vert}
+ E_{\rm xc}[\rho]$  

그리고, 위 식을 KS방정식의 고유값($\varepsilon_i$)과 전자밀도로 나타낼 수 있다. 즉, 식 (6)에 $\psi_i^{\ast}$ 을 곱하고 적분을 한 후, $i$에 대하여 합을 취하면 식 (10) 은 다음의 식으로 바뀐다.
 

\begin{displaymath}
E[\rho] = \sum_i^N \varepsilon_i - \frac12 \int d {\bf r}\in...
...ilon_{\rm xc}(\rho({\bf r}))- V_{\rm xc}(\rho({\bf r}))\right]
\end{displaymath} (11)


 

그러나, 교환-상관에너지 범함수의 정확한 형태를 알지 못하기 때문에, 그 형태를 적당히 근사를 하여야 한다. Kohn과 Sham은 국소밀도근사법(LDA)을 제안하였다. 국소밀도근사는 간단하면서도 널리 사용되는 근사로서, 어떤 한 점에서 전자 하나가 가지는 교환-상관 에너지를 그 점에서 같은 밀도를 가지는 균일한 전자 기체에서 하나의 전자가 가지는 교환-상관 에너지와 같다고 생각을 하는 것이다. 즉, 교환-상관 에너지가 국소적인 성격을 가진다고 가정한다. 이러한 단순한 가정에도 불구하고 그동안 LDA는 상당히 성공적으로 원자, 분자, 고체의 성질을 계산하였다.

하지만, 근본적으로 LDA는 주위 밀도의 변화를 무시되는 근사이므로, 분자나 고체 표면등과 같이 전자 밀도가 심하게 변하는 계에서는 잘 맞지 않을 것을 예상할 수 있다. 실제로 계산을 해 보면 다음과 같은 경향이 있음이 알려져 있다. 즉, 실험과 비교하여 분자나 고체의 결합 에너지를 크게 기술하며, 약하게 결합된 분자와 고체의 결합 길이와 결정 상수를 짧게 기술한다는 것이다[3]. 이런 LDA의 한계를 극복하기 위하여 교환-상관 에너지를 계산할때 밀도의 국소값뿐 아니라, 그것의 미분값도 반영하는 근사를 생각하게 되었는데, 이를 Generalized Gradient Approximation(GGA) 라고 부른다.

 

2. 슈도퍼텐셜

원자는 속전자와 가전자로 이루어져 있다. 그런데, 원자들이 결합을 하여 고체를 이룰때, 속전자는 원자핵근처에 강하게 구속되어 있어, 가전자만이 실제 고체결합에 관계된 것으로 생각할 수 있다. 그래서, 고체를 원자핵과 속전자가 주기적으로 배열되고 있고, 그 주위로 구속되어 있지 않은 가전자들이 움직이는 것으로 생각한다. 이와 같은 모형에서는 가전자끼리의 상호 작용과, 가전자와 원자핵과 속전자를 합한 이온과의 상호 작용이 있다. 첫번째 경우는 하트리나 DFT로 근사를 하며, 두번째는 슈도퍼텐셜 방법을 이용한다.

슈도퍼텐셜은 가전자의 성질을 효과적으로 기술하면서, 강한 전자-이온 퍼텐셜을 부드럽게 바꾼것이다. 그렇게 함으로써, 슈뢰딩거방정식을 적은 수의 평면파로 전개하여 풀 수 있게 된다. 이와 비슷한 접근은 1930년대에 Fermi에 의해서 시도되었다. 그 후 Herring[4]가 Orthogonalized plane-wave 방법(OPW)을 그리고, Phillips 과 Kleinman [5]이 OPW 의 직교항은 척력 퍼텐셜을 준다는 것을 보였다. 그 후 실험값으로부터 슈도퍼텐셜을 얻어서 여러가지 계산을(EPM 방법) 하였다[6]. 1980년대에 들어와서, Hamann, Schlüter, Chiang[7]은 제일 원리적, 즉 실험값을 사용하지 않고 원자번호만을 입력하여 슈도퍼텐셜을 만들어 내는 방법을 제안하였다.

대부분의 전자구조 계산에서 사용되는 슈도퍼텐셜은 모든 전자를 고려한 원자계산의 결과를 이용하여 만든다. DFT를 바탕으로 한 슈도퍼텐셜은 구면 가리기 근사를 사용하여, 다음의 반경 KS 방정식을 자체 총족적으로 풀어서 구한다.
 

\begin{displaymath}
\left[ \frac{-1}{2} \frac{d^2}{d r^2} + \frac{l(l+1)}{2 r^2}...
...rho;r] \right] r R_{nl}(r) = \varepsilon_{nl}   r R_{nl}(r),
\end{displaymath} (12)


 

여기서 $V[\rho;r]$은 자체 충족적으로 구한 전자 하나가 느끼는 퍼텐셜로 다음과 같다.
 

\begin{displaymath}
V[\rho;r] = \frac{-Z}{r}+ V_{\rm H}[\rho;r] + V_{\rm xc}[\rho(r)],
\end{displaymath} (13)


 

여기서, $\rho(r)$은 채워진 파동함수$R_{nl}(r)$의 전자 밀도 합이며, $V_{\rm H}[\rho;r]$은 하트리 퍼텐셜을, $V_{\rm xc}[(\rho(r))]$는 교환-상관 퍼텐셜을 나타낸다. 대부분 슈도퍼텐셜은 다음의 4가지 조건을 만족하도록 만들어진다. 첫째는, 슈도퍼텐셜로 부터 구한 슈도파동함수는 마디를 가지지 않는다. 두번째로는 규격화된 슈도파동함수(PP)는 우리가 정한 어떤 $r_c$바깥에서는 모든 전자를 고려한 파동함수(AE)와 일치한다.
 

\begin{displaymath}
R_l^{\rm PP}(r) = R_l^{\rm AE}     \mbox{for\hspace{0.2cm}}  r> r_c \
\end{displaymath} (14)


 

세번째로, $r_c$안에 있는 두 파동함수의 전하량은 같다.
 

\begin{displaymath}
\int_0^{r_c}\vert R_l^{\rm PP}\vert^2 r^2 d r = \int_0^{r_c}\vert R_l^{\rm AE}\vert^2 r^2 d r
\end{displaymath} (15)


 

그리고, 끝으로 실제와 슈도원자의 고유값이 일치한다.
 

\begin{displaymath}
\varepsilon_l^{\rm PP} = \varepsilon_l^{\rm AE}
\end{displaymath} (16)


 

만약 슈도퍼텐셜이 위의 4가지을 만족한다면, 우리는 이것을 ``크기를 보존하는 슈도퍼텐셜''이라고 부른다[7]. 위의 조건을 만족시키면서 슈도퍼텐셜을 만드는 방법은 여러가지가 있는데 [7,8,9,10], 우리가 사용하는 것은 Troullier 와 Martins의 방법이다[10].

최근에는 위의 조건에서 세번째 조건을 제거한, 즉 슈도파동함수의 크기에 대한 구속조건을 완화한 초연성 슈도퍼텐셜방법도[11] 많이 사용되고 있다. 이것은 첫번째 줄에 있는 원자나 3$d$ 전자를 포함하는 전이 금속과 같이 매우 국재화된 전자를 가지고 있는 계에 주로 적용된다.

 

3. 총에너지 공식

Korn-Sham 방정식을 푸는 방법으로 여러가지가 있지만, 슈도퍼텐셜을 사용할 경우 특히 유용한 방법은 전자 파동함수를 역공간에서 평면파를 기저로 전개를 하여 푸는 것이다. 평면파를 사용하면, 몇가지 이점이 있다. 우선 해밀토니안 행렬을 쉽게 계산 할 수 있다. 그리고, 평면파의 갯수를 늘려서 우리가 원하는 정확도까지 얻을 수 있다. 특히, 총에너지와 Hellman-Feynman 힘이 평면파에서는 간단한 형태로 표현이 된다[12].

Bloch 정리에 따르면, 파동함수는 다음과 같이 쓰여진다.
 

\begin{displaymath}
\psi_{n {\bf k}}({\bf r}) = \frac{1}{\sqrt{V}} \sum_{\bf G}
\psi_{n{\bf k}}({\bf G}) e^{i({\bf k+G})\cdot {\bf r}} ,
\end{displaymath} (17)


 

여기서, $V$는 전체 결정의 부피, $n$은 밴드 표시첨자, k는 브릴루앙 영역의 한 점을 나타낸다. 그리고, 다른 양들을 평면파로 전개를 하면 다음과 같다.
 

$\displaystyle \rho({\bf r})$ $\textstyle =$ $\displaystyle \sum_{\bf G} \rho({\bf G}) e^{i{\bf G}\cdot {\bf r}},$ (18)
$\displaystyle V_H({\bf r})$ $\textstyle =$ $\displaystyle \sum_{\bf G} V_H({\bf G}) e^{i{\bf G}\cdot {\bf r}},$ (19)
$\displaystyle V_{\rm xc}({\bf r})$ $\textstyle =$ $\displaystyle \sum_{\bf G} V_{\rm xc}({\bf G}) e^{i{\bf G}\cdot {\bf r}},$ (20)


위에서
G 는 역격자 벡터이고, $\rho({\bf G})$, $V_H({\bf G})$, $V_{\rm xc}$ 는 각각 전자밀도, 하트리 에너지, 교환 상관 퍼텐셜의 푸리에 변환의 성분이 된다. Kohn-Sham방정식(식 6)에 위의 식들은 대입한 후 정리를 하면, 다음과 같은 역격자 공간에서 단일 전자 방정식을 구할 수 있다.
 

\begin{displaymath}
\sum_{{\bf G'}}\left[ \frac12({\bf k}+{\bf G'})^2 \delta_{{\...
...({\bf G'})
= \varepsilon_{n {\bf k}}\psi_{n {\bf k}}({\bf G}).
\end{displaymath} (21)


 

그리고, 식(10)에 이온-이온 에너지를 포함하여, 총에너지를 다음과 같은 푸리에 계수로 나타낼 수 있다.
 

$\displaystyle E_{\rm tot}$ $\textstyle =$ $\displaystyle T+E_{\rm ion}+E_{\rm H}+E_{\rm xc}+ E_{\rm ion-ion}$ (22)
  $\textstyle =$ $\displaystyle \left[\frac12 \sum_{n {\bf k G}} \vert\psi_{n {\bf k}}({\bf G})\vert^2
({\bf k}+{\bf G})^2 \right]$  
  $\textstyle +$ $\displaystyle \left[ \Omega \sum_{{\bf G}} V^{\rm L}({\bf G})\rho({\bf G})
+ \s...
...{\bf G}) \psi_{n {\bf k}}({\bf G'})
V_l^{\rm NL}({\bf k+ G},{\bf k+G'}) \right]$  
  $\textstyle +$ $\displaystyle \left[\frac{\Omega}{2} \sum_{{\bf G}}V_{\rm H}(-{\bf G}) \rho({\b...
...)\right]
+\left[ \sum_{{\bf G}} \epsilon_{\rm xc}({\bf G}) \rho({\bf G})\right]$  
  $\textstyle +$ $\displaystyle \left[ \frac12 \sum_{\bf R \ne R'} \frac{Z Z'}{\vert{\bf R - R'}\vert}\right] .$  
       

여기에서, $V^{\rm L}$은 국소 슈도퍼텐셜을 $V_l^{\rm NL}$은 비국소 슈도퍼텐셜을 나타낸다. 그리고, $\Omega$는 결정의 총부피를, $\epsilon_{\rm xc}({\bf G})$는 교환-상관에너지의 푸리에 계수를 나타낸다. 그러나, 위의 식(22)에는 $E_{\rm ion-ion}$, $V_{\rm H}({\bf G}=0)$, 그리고, $V^{\rm L}({\bf G}=0)$과 같은 발산하는 양이 있다. 하지만, 이 양은 서로 상쇄됨을 보일 수 있다. 그래서 다음과 같이 바뀐다.
 

\begin{displaymath}
E_{\rm tot} = T + E_{\rm ion}^{'} + E_{\rm H}^{'} +E_{\rm xc}+
N(Z \alpha_1 + \gamma_{\rm Ewald}).
\end{displaymath} (23)


 

여기서, 프라임기호($'$)가 붙은 항은 합을 할 때 발산하는 항을 빼고 계산을 한다는 것이다. 그리고, $\alpha_1$$\gamma_{\rm Ewald}$는 다음과 같이 정의된다.
 

\begin{displaymath}
\alpha_1 = \lim_{{\bf G}\rightarrow 0} \left[ V_{\rm ion}^{\...
...t[V_{\rm ion}^{\rm L}({\bf r})+
\frac{Z}{r}\right] d {\bf r},
\end{displaymath} (24)


 

\begin{displaymath}
\gamma_{\rm Ewald} = \frac12\left[\sum_{\mu}^{'} \frac{Z^2}{...
...rt}
- \frac{1}{\Omega_{at}}\int \frac{Z^2}{r} d {\bf r}\right]
\end{displaymath} (25)


 

여기서 $\Omega_{at} = \Omega/N$이다. 그리고, 이것을 Kohn-Sham 고유값이 들어간 식으로 바꾸어 다음과 같은 최종형태를 얻는다.
 

$\displaystyle E_{\rm tot} (\mbox{per atom})$ $\textstyle =$ $\displaystyle \sum_{n \bf k}\varepsilon_{n \bf k}
-\frac{\Omega_{at}}{2} \sum_{{\bf G}\neq 0} V_{\rm H}(-{\bf G})\rho({\bf G})$ (26)
  $\textstyle +$ $\displaystyle \Omega_{at} \sum_{{\bf G}} \left[ \epsilon_{\rm xc}({\bf G})
- V_...
... xc}({\bf G})\right] \rho({{\bf G}})
+ Z_{\rm c} \alpha_1 + \gamma_{\rm Ewald}.$  

 

참고 문헌

 
1
P. Hohenberg and W. Kohn, Phys. Rev.136, B864 (1964)
2
W. Kohn and L. J. Sham, Phys. Rev. 140, A1113 (1965)
3
R. O. Jones and O. Gunnarson, Rev. Mod. Phys. 61, 689(1989)
4
C. Herring, Phys. Rev. B 57, 1167 (1940).
5
J. C. Phillips and L. K. Kleinman, Phys. Rev. B 116, 287 (1959)
6
M. L. Cohen, Phy. Rep. 110, 293 (1984)
7
D. R. Hamann, M. Schlüter, and C. Chiang, Phys. Rev. Lett. 43, 1494 (1979)
8
G. P. Kerker, J. Phys. C 13 L189 (1980)
9
G. B. Bachelet, D. R. Hamann, and M. Schlüter, Phys. Rev. B 26, 4199 (1982)
10
N. Troullier and J. L. Martins Phys. Rev. B 43, 1993 (1991)
11
D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
12
J Ihm, Alex Zunger, and M. L. Cohen, 1979, J. Phys. C 12, 4409; J Ihm, Rep. Prog. Phys. 51 (1988)