VASP (1) Brillouin Zone Sampling

VASP 계산의 이론적 배경에 대해 공부한 내용을 작성하려고 한다. 이 내용은 VASP 개발자 중 한 명이 진행한 workshop 강의를 토대로 한다. 내용이 어려워 제대로 이해하지 못한 부분들이 많지만, 최대한 내가 이해한 바를 쓰려고 한다.

https://youtu.be/vJkNv095Aj8

https://youtu.be/dVguc4W1hrQ

 

아래는 비교적 최근에 진행된 비대면 강의이다.

https://youtu.be/Fv3F4LHGPuc

https://youtu.be/v7gc98lG6Wo

 

사실 이 내용들은 VASP에만 해당하는 건 아닌 것 같고, Quantum Espresso 등 다른 DFT 계산 프로그램에서도 사용되는 원리인 것 같다.

 

1. Sampling Brillouin zone

Periodic boundary condition을 갖는 결정(crystal)의 물리량들은 1st Brillouin zone(BZ)에서의 적분으로 계산될 수 있다.

$$\begin{equation}\tag{1}f(r)=\frac{1}{\Omega_{BZ}}\int F(\mathbf{k})d\mathbf{k}\end{equation}$$

여기서 $\Omega_{BZ}$는 BZ의 부피이다. 예를 들면 charge density는 다음과 같이 계산된다.

$$\rho(\mathbf{r})=\frac{1}{\Omega_{BZ}}\sum_n \int_{BZ}f_{n\mathbf{k}}|\psi_{n\mathbf{k}}(\mathbf{r})|^2 d\mathbf{k}$$

이때 $f_{n\mathbf{k}}$는 occupancy에 관한 식이다. 적분으로 표현된 1번 식을 다음과 같이 discrete k-point에 대한 weighted sum으로 바꾸어 계산할 수 있다.

$$\begin{equation}\tag{2}f(\mathbf{r})=\sum_{\mathbf{k}}w_\mathbf{k}F(\mathbf{k})\end{equation}$$

이렇게 하기 위해서는 BZ 안에서 k-point를 골라야 한다. 이것을 BZ sampling, 또는 k-point sampling이라 한다. 당연하게도 k-point를 더 많이 고를수록 더 정확한 계산을 할 수 있다. 하지만 많은 k-point를 고르면 계산 비용이 커지기 때문에, k-point 수를 늘려 가며 계산한 결과값이 어느정도 수렴하는 지점이 있다면 그 지점의 k-point 개수를 쓰는 것이 좋다. 이때 BZ 내에서 k-point를 균일(uniform)하게 고르는 것이 빠른 k-point convergence에 좋다고 한다. BZ sampling의 간단한 예시로, 2D cubic lattice에서 $4\times4$ mesh를 사용할 경우 다음과 같이 k-point를 고를 수 있다.

즉 $4\times4$라는 의미는 각 reciprocal lattice vector를 따라 4개의 k-point를 고른 것을 의미한다. 3D lattice에서는 $N_1\times N_2\times N_3$로 k-point mesh를 표현한다. 그리고 이와 같이 균일한 BZ sampling을 하기 위해 널리 사용되는 것이 Monkhorst-Pack(MP) grid(또는 mesh)이다. 그런데 BZ를 잘 보면, 표시된 삼각형 부분을 돌리고 뒤집고 하면(symmetry operation) 모든 BZ에 딱 맞게 들어간다는 것을 알 수 있다. 이렇게 해당 lattice가 가지는 symmetry operation을 해서 전체 BZ를 나타낼 수 있는 가장 작은 영역을 irreducible Brillouin zone(IBZ)라 한다. 그렇다면 이 IBZ 안의 k-point에 대해서만 계산한 뒤, 각 k-point와 동등한 점들이 BZ 전체에 몇 개가 있는지를 세서 적절한 weight를 곱하고 더해주면 BZ 전체의 k-point에 대해 sum을 한 셈이나 다름없다. 위 그림을 예로 들면 전체 BZ 안에서 $k_1$와 동등한 점은 자기 자신을 포함해서 총 4개, $k_2$는 4개, $k_3$은 총 8개가 있다. 그러면 BZ 안의 전체 k-point 수가 16개이므로 각각의 weight는 4/16, 4/16, 8/16이 되고, 따라서 2번 식은 다음처럼 된다.

$$\frac{1}{4}F(\mathbf{k_1})+\frac{1}{4}F(\mathbf{k_2})+\frac{1}{2}F(\mathbf{k_3})$$

K-point mesh를 만드는 알고리즘은 일반적으로 다음과 같다고 한다.

 

• 동일한 간격의 mesh를 만든다.
• 필요하다면 mesh를 옮긴다(KPOINTS 파일에서 설정할 수 있다). 
• Bravais lattice에 따른 symmetry operation을 모든 k-points에 대해 적용한다.
• IBZ를 결정한다.
• 적절한 weight를 계산한다.

K-point의 개수는 cell의 dimension과 밀접한 관련이 있다. Real lattice의 lattice parameter가 클수록 reciprocal lattice의 lattice vector 길이는 줄어들기 때문에, BZ 또한 작아지게 된다. 따라서 작은 cell일 때보다 더 적은 개수의 k-point로도 충분하다. 극단적인 예로 무한한 cell이라면 reciprocal space에는 원점만 있을 것이므로 그 한 점에 대해서만 계산해도 되는 것이다. 또한 한 방향으로 길쭉한 cell인 경우, 예를 들어 $z$축 방향으로의 lattice constant가 다른 두 축보다 훨씬 크다면, 균일한 BZ sampling을 하기 위해서는 $z$축에 해당하는 reciprocal lattice vector를 따라서 sampling해야 하는 k-point 개수가 더 줄어든다.

 

어떤 축으로의 k-point 개수(subdivision)는 그 축 방향으로의 periodic boundary condition을 얼마나 고려하고 싶은지를 결정하는 것이다. 따라서 큰 subdivision을 사용하면 해당 축 방향으로의 interaction을 더욱 정확히 기술하게 된다. 모든 축 방향으로 무한한 주기성을 갖는 bulk crystal의 경우에는 앞서 언급한 cell의 dimension과 k-point convergence를 확인하여 적절한 subdivision을 사용해야 한다. 하지만 분자 하나의 에너지를 계산하고 싶을 때와 같이 periodic boundary condition의 영향을 최소하하고 싶을 때도 있다. 이럴 때는 분자 주변에 vacuum을 충분히 설정하여 cell의 크기를 충분히 크게 한 다음 k-point mesh를 $1\times 1\times 1\times$로 하여 분자 사이의 interaction을 최소화해야 한다.

분자에 periodic boundary condition가 적용되면 이렇게 분자들이 periodic하게 배열되기 때문에 분자 주변에 충분히 vacuum을 설정한 다음 $1\times1\times1$의 mesh를 사용해야 분자 사이의 interaction을 최소화할 수 있다.

또한 slab을 만들어 surface 특성을 계산하고 싶을 때는 slab의 수직한 방향(예를 들어 $z$축)으로 vacuum을 충분히 설정한 다음 $z$축에 해당하는 subdivision을 1로 해야 한다. 즉 $N_1\times N_2\times 1$의 mesh를 사용해야 slab과 slab 사이의 interaction을 최소화할 수 있다.

Slab에 periodic boundary condition이 적용되면 이와 같이 slab이 쌓여 있는 형태가 되기 때문에 slab과 slab 사이의 vacuum을 충분히 만들어 주고 해당 subdivision을 1로 해야 slab 사이의 interaction을 최소화할 수 있다.

2. $\Gamma$-centered vs. around-$\Gamma$

한편 k-point mesh는 reciprocal lattice의 원점($\Gamma$ point)를 포함하는 경우($\Gamma$-centered)와 아닌 경우(around-$\Gamma$)로 나눌 수 있다. 

MP mesh를 만들 때 짝수 개의 k-point를 지정하면 경우 왼쪽 그림처럼 around-$\Gamma$ mesh가 만들어지고, 홀수 개의 k-point를 지정하면 오른쪽 그림처럼 $\Gamma$-centered mesh가 만들어진다. 위 경우에는 두 가지 모두 lattice의 symmetry를 깨뜨리지 않기 때문에 어떤 mesh를 사용하더라도 큰 문제는 없지만, 그렇지 않은 경우도 있다. 예를 들어 hexagonal lattice를 살펴보자.

Hexagonal lattice에서 왼쪽과 같이 around-$\Gamma$ mesh를 만들면, 이 mesh는 hexagonal lattice가 가지고 있는 60도 회전 대칭을 만족시키지 못한다. 그런데 위에서 설명한 것처럼 mesh를 만드는 과정에서 lattice의 symmetry operation을 하는데, 그러면 가운데 그림처럼 새로운 k-point들이 생기면서 원래 의도했던 것과는 다른 mesh가 만들어진다. 왜냐하면 원래 mesh가 60도 회전 대칭을 만족시키지 못하기 때문에 k-point들이 서로 겹쳐지지 않기 때문이다. 그 결과 mesh가 균일하지 못하게 되고, 그러면 k-point convergence가 잘 되지 않는다. 이를 방지하기 위해서는 오른쪽 그림처럼 $\Gamma$-centered mesh를 구성해야 한다. VASP wiki에 따르면 face-centered cubic(FCC), hexagonal, 그리고 face-centered orthorhombic lattice system에서는 $\Gamma$-centered mesh를 써야 한다. 

 

위에서 언급한 대칭성이 깨지는 경우가 아니더라도, $\Gamma$를 포함하느냐 아니냐는 IBZ 안의 k-point 개수에 영향을 준다. 예를 들어 cubic lattice의 경우, $7\times7\times7$ MP mesh를 만들면 IBZ 안에 20개의 k-point가 선택되는데, $8\times8\times8$ MP mesh를 만들어도 20개가 선택된다. 그래서 MP mesh로 할 경우 subdivision을 짝수($\Gamma$ 미포함)로 할 때와 홀수($\Gamma$ 포함)로 할때의 k-point convergence를 따로 확인하는 것이 좋다고 한다. 왜냐하면 k-point convergence를 결정하는 것은 IBZ 안에서 실제로 선택된 k-point의 개수인데, subdivision의 크기를 늘리는 것이 항상 IBZ 안의 k-point 개수를 증가시키는 것이 아니기 때문이다. 사실 지금껏 언급한 여러 문제들을 신경쓰지 않고 싶으면 항상 $\Gamma$-centered mesh를 쓰면 되긴 한다. VASP 강의에서도 강의자 본인은 항상 $\Gamma$를 포함시켜서 계산한다고 말했다. 하지만 system에 따라, 또 계산하고자 하는 물성에 따라 k-point convergence를 항상 잘 확인할 필요가 있다.

 

3. Smearing

Metal의 경우, 0K에서는 BZ 내에서 전자의 occupancy가 불연속적으로 변하는 경계가 있다. 왜냐하면 metal은 valence band가 꽉 채워져 있지 않기 때문인데, 이러한 경계면을 Fermi surface라 한다. 이렇게 불연속적인 함수가 있을 때는 1번 적분식을 2번 sum으로 근사하려면 많은 k-point가 필요하고, k-point convergence가 어려워진다. 이를 해결하기 위해 smearing이라는 방법을 쓰는데, Fermi suface를 흐릿하게 만들어 불연속적인 경계를 부드럽게 만드는 것이다. 하지만 너무 많이 smear하게 되면 당연히 부정확한 계산 결과를 얻을 수 있기 때문에, 적절한 smear parameter를 선택해야 하고 그에 따른 k-point convergence 또한 확인해야 한다. 그리고 ISMEAR = -4 or -5에서는 항상 $\Gamma$-centered mesh를 사용해야 한다고 한다.

 

추가 참고자료:

Density functional theory: A practical introduction

https://www.vasp.at/wiki/index.php/KPOINTS

KPOINTS - Vaspwiki

https://sites.psu.edu/dftap/2019/02/01/shifting-of-kpoints-in-hexagonal-lattices/

Shifting of k points in hexagonal lattices | Density Functional Theory and Practice Course