VASP 계산 경험적 지식 및 유용한 사이트 모음

VASP 계산을 하면서 경험적으로 얻은 지식과 참고하기 좋은 사이트들을 정리하려고 한다. 아직 경험이 많은 게 아니기 때문에 지속적으로 업데이트하려고 한다. 

 

1. Electronic SCF calculations

Electronic scf step의 알고리즘은 ALGO tag로 정해주게 되는데, 보통 많이 사용되는 방법이 Fast, Normal, All 정도가 있다. 이 외에 하이브리드 functional (e.g. HSE06) 의 경우 Damped를 사용하기도 한다 (사실 이 경우에는 Fast나 Normal을 사용할 수 없다). 각 방법들에 대해 자세히 설명하기에는 너무 길어서 여기에 적지 않고 [1]을 참고하길 바란다. Practical한 측면에서 말하자면 속도는 일반적으로 Fast > Normal > All 순이며 robustness는 반대 순서로 생각하면 된다. 그래서인지 VASP default ALGO는 Normal이다. 한편 Magnetic moment를 고려하게 되는 spin-polarized 계산에서, 특히 GGA+U에서는 Fast나 Normal이 가끔 잘 수렴하지 않는 경우가 있어서, 개인적으로는 GGA+U를 사용하는 경우나, magnetic moment 값이 크거나 중요한 system의 경우에는 All을 사용하고 그렇지 않을 때는 Fast를 사용하는 것이 더 편했다.

 

Fast나 Normal을 사용할 때 조금 더 보수적이지만 안정적으로 진행하고 싶다면 mixing parameter를 조절해주면 좋다. Vaspwiki 등 공식 자료에서도 Fast나 Normal의 수렴이 어려운 경우 (예를들면 slabs, magnetic systems, molecules, clusters) AMIX=0.2, BMIX=0.001 같은 값을 사용해보라고 하고 있다. 

 

2. 구조최적화

구조최적화 알고리즘은 IBRION으로 정해주게 되고, 1, 2, 3의 3가지가 있다. 역시 자세한 알고리즘의 설명은 여기서 하지 않고 [1]을 참고하길 바란다. 내가 이해한 바를 간략히 설명하자면 1과 2의 formalism은 potential energy surface (PES)가 quadratic 한 상태, 즉 어느정도 local minimum 근처에 있을 때를 가정하고 있다. 그래서 초기 구조가 local minimum 근처에 있거나 equilibrium structure와 비슷하다면 수렴 속도가 빠르다. 반면 초기 구조가 local minimum과 매우 멀리 떨어져 있는 경우라면 불안정한 알고리즘이다. 그렇지만 1에 비해 2는 quadratic이 아니라도 local minimum으로 어느정도 잘 가는 편이라고 한다. 3은 초기 구조가 local minimum과 멀리 있을 때 가장 사용할 만한 알고리즘이지만, 속도는 느린 편이다. VASP wiki에는 3을 잘 사용하면 2보다 빠르다고 하는데, 후에 설명하듯이 그렇게 하려면 손이 많이 가고 경험이 많이 필요할 것 같다. 일반적으로는 속도는 1 > 2 > 3 이고 안정성은 3 > 2 > 1 정도로 볼 수 있겠다. 한편 1과 3은 POTIM을 잘 설정해줘야 할 필요가 있고, 특히 1은 NFREE를 설정해주지 않으면 나의 경우 memory 부족 에러가 떠서 이것까지 설정하려면 손이 많이 간다. 3은 SMASS 값을 설정해줄 수도 있는데, 이것 또한 수렴성에 영향을 준다. 따라서 속도와 안정성, 그리고 parameter 설정의 편리성을 모두 고려하면 2를 쓰는 것이 무난해 보인다. 그래서인지 일반적으로 2가 많이 사용되는 편이다. 하지만 앞에서 말한대로 초기 구조의 상태가 나쁘거나, 어느정도 수렴한 시스템을 더욱 tight하게 수렴시키는 등의 목적이라면 3 및 1을 적절히 활용해볼 수 있고, 이에 대한 예시는 [2]에 있다. 근데 위에서 말했듯이 1이나 3은 손이 많이 가고, 또 3은 표면 계산에서는 vacuum 때문인지 표면이 무너지는 일이 발생하는 등 이러나 저러나 2가 제일 무난하게 쓸만한 것 같다. 하지만 1과 3이 거의 안쓰이는 건 아니고, 나중에 언급할 NEB 계산에서 주로 쓰이게 된다. 미리 말하자면 NEB 계산에서는 2를 쓰기 어렵다. 

 

2.1. Preconverge with low precision → reconverge with high precision

구조최적화 시 초기 구조는 아무래도 equilibrium 구조와는 좀 떨어져 있을 수밖에 없고, 이럴 때 높은 precision의 electronic scf 계산을 굳이 할 필요가 없는 경우가 많다. 이 부분은 Henkelman group의 포럼[4], 그리고 [5]에서도 언급되고 있는데, 예를 들면 ENCUT이나 KPOINTS를 처음부터 높게 설정해서 계산할 필요가 없다는 뜻이다. 대신 낮은 ENCUT (나같은 경우는 ENMAX로 맞췄다)과 낮은 밀도의 KPOINTS grid (예를 들면 Gamma point) 를 설정해서 preconverge시키고, 이 구조로부터 다시 ENCUT과 KPOINTS를 늘려서 계산을 하는 것이 더 효율적이다. 이렇게 하면 구조최적화를 두 번 하는거라 비슷하거나 더 느릴 수도 있다고 생각할 수도 있겠지만, 처음부터 높은 ENCUT과 KPOINTS로 구조최적화할 때보다 경험상 적어도 2배 이상 빠르게 수렴하였다. 그래서 난 기본적으로 이런 방식으로 구조최적화를 진행하고 있다.

 

2.2. Symmetry constraint

INCAR에서 ISYM tag는 초기 구조가 가지고 있는 symmetry를 깨뜨리지 않게끔 구조에 symmetry constraint을 가할지 말지를 결정한다. Symmetry를 이용하면 계산 속도 측면에서는 장점이 있지만 종종 symmetry constraint으로 인해 더 안정한 구조로 수렴하지 못하는 경우가 있다. 나는 개인적으로 lattice parameter 최적화에는 ISYM=2로 계산하고 (주로 데이터베이스에서 가져온 구조파일에서 시작라는 경우) 그 이후 defect나 adsorption 계산 등에서는 ISYM=0으로 계산하는 편이다.  

 

2.3. Pulay stress

셀의 부피도 변하는 구조최적화 (e.g. ISIF=3)의 경우, ENCUT으로 basis set 크기를 지정하더라도 셀의 부피가 변하기 때문에 역격자 공간의 크기가 변하게 되고, 이로 인해 실질적인 ENCUT이 달라지는 효과가 생긴다. 이로 인해 내가 설정한 ENCUT에서 완전히 수렴한 경우보다 부피를 작게 만드는 경향이 생긴다고 하는데, 이를 Pulay stress라고 한다. 실제로 ISIF=3의 계산에서, 한번 수렴시킨 구조를 동일한 input files를 사용해서 또 수렴시키면, 이미 수렴된 구조라서 ionic step 한번으로 끝나는 게 아니라 ionic step을 더 밟는 경우가 많이 있다. 그래서 내가 설정한 ENCUT 및 수렴 조건 안에서 제대로 수렴시키기 위해 구조최적화를 여러 번 반복할 수도 있다. 특히 lattice parameter가 중요한 상황이라면 이렇게 해볼 수 있겠다. 구조최적화를 그렇게까지 엄밀하게 하지 않더라도, 부피가 변한 경우에는 실질적인 ENCUT이 달라졌기 때문에 최적화된 구조의 total energy를 정확히 구하려면 single-point 계산을 한번 더 해야 한다. 그렇게 해야 내가 설정한 ENCUT에서의 에너지를 얻을 수 있다. 여기에 대한 자세한 내용은 Energy vs volume Volume relaxations and Pulay stress - VASP Wiki 및 [1]을 참고하면 된다.

 

3. Frequency 계산

Frequency 계산은 NSW=1; IBRION=5, POTIM=0.01을 사용하고 있다. 여기에 대해서는 특별히 공부를 해 보거나 계산을 테스트해본 적은 거의 없다. Vaspwiki에는 frequency 계산을 할 때 IBRION을 5 또는 6으로 할 수 있다고 나오는데, 6을 할 경우 selective dynamics를 할 수 없다. 또 NFREE라는 tag를 설정해줄 수도 있는데, vaspwiki에는 설명이 제대로 안 나와 있는 것 같은데 OUTCAR를 직접 살펴봤을 때 IBRION=5인 경우 NFREE의 default가 2인 것으로 확인했다. 주로 이 default를 사용하는 것 같다. 또한 frequency 계산을 할 때에는 force를 정확히 구하는 것, 즉 EDIFF를 tight한 조건으로 수렴시키는 것이 좋다고 하는데, 엄밀한 계산을 요하지 않는 경우라면 그냥 구조최적화할 때와 같은 조건으로 계산해도 큰 문제는 없는 것 같다. Phonon계산을 하는 것 같으면 신경써야 할 것들이 많아 보이는데 내 연구실에서는 phonon 계산을 하지 않아서 크게 어려운 부분은 없는 것 같다.

 

한편 IBRION=5 에서는 NCORE=1로 설정하는 것이 일반적이나, 내 경우 이렇게 하면 segmentation fault error가 떴다. 메모리 부족으로 인한 오류인 것 같아 NCORE를 늘려서 계산하려고 하니 이번엔 VASP 내부적인 에러가 떴다. 이는 https://mattermodeling.stackexchange.com/questions/9004/ncore-bigger-than-1-when-performing-phonon-vibrational-calculationibrion-5-in 를 참고해서 해결할 수 있었다. 즉 ISYM=0으로 놓고 하면 해결된다.

 

4. NEB 계산

NEB 계산은 transition state (TS)를 찾는 방법 중 하나이다. NEB 계산의 알고리즘은 다음 영상에서 잘 설명되어 있다.

https://youtu.be/2PotR4_sV_0

 

아래에 설명할 NEB 계산에 대한 내용의 대부분의 출처는 Graeme Henkelman group의 forum ([4])이다. NEB 계산은 initial state (IS)와 final state (FS)를 각각 구조최적화한 뒤, images로 이어진 band를 relax 하는 과정으로 이루어진다. IS와 FS의 구조최적화는 특별한 건 없고 1번에서 설명한 것과 같은 일반적인 방법들을 따르면 된다. 물론 이 구조최적화와 band relaxation은 같은 precision (ENCUT, KPOINTS 등)에서 이루어져야 한다. NEB 계산에서는 ISYM=0을 하는 것이 권장된다. Band를 잘 수렴시키고 TS를 잘 찾기 위해서는 먼저 IS와 FS를 local minimum에 잘 위치시키는 것이 중요한데, symmetry constraint으로 인해 이것이 제한될 수도 있다. 또한 이후 진행되는 band relaxation에서도 ISYM=0을 해서 symmetry constraint으로 인해 band가 제대로 수렴되지 않는 것을 방지하도록 한다. 그리고 NEB 계산에서도 낮은 precision에서 수렴시킨 구조로부터 시작해서 높은 precision으로 다시 수렴시키는 방식이 유효하고, 전체적인 시간 역시 단축된다. 즉 IS 및 TS 구조최적화와 band relexation을 모두 낮은 precision으로 진행한 다음, 이렇게 해서 얻은 band에서 높은 precision으로 다시 시작하는 것이다. 물론 높은 precision으로 재수렴시킬 때는 IS와 FS 역시 높은 precision으로 수렴시킨 구조여야 한다. 

 

4.1. Climbing-image NEB

Climbing-image NEB (CI-NEB)는 기존의 NEB 알고리즘을 조금 수정하여 이미지들 중 가장 높은 에너지를 갖는 이미지를 saddle point로 수렴시키게끔 하는 방법이다. 이에 대한 설명 역시 위 영상에서 볼 수 있다. 기존의 NEB는 TS를 잘 guess를 위한 방법이라고 한다면, CI-NEB에서는 가장 높은 에너지를 갖는 이미지 그 자체가 TS가 되게끔 하는 것이다. 모든 이미지가 다 원하는 force criterion (EDIFFG)에 수렴하지 않더라도 가장 에너지가 높은 이미지가 수렴하면 그걸 TS로 생각해도 큰 문제는 없는 것 같다. (아니면 Dimer method 또는 improved dimer method를 사용해서 TS를 더 자세하게 찾아볼 수도 있다.) 이러한 장점 때문에 거의 대부분 CI-NEB를 사용한다. 

 

4.2. Optimizer 선택

NEB 계산은 optimizer 선택이 꽤 중요할 수 있다. 보통 구조최적화를 할 때는 IBRION=2를 많이 선택하지만, NEB 계산에서는 이것을 쓸 수 없다. 이걸 써도 에러가 뜬다거나 계산이 무조건 터진다는 건 아니지만, 수렴성이 나쁠 가능성이 크다. 왜냐하면 NEB 계산은 NEB force를 relax하는 방식인데, IBRION=2는 에너지를 낮추는 방식으로 작동하기 때문이다. 그래서 force 기반 optimizer인 IBRION=1 또는 3을 사용해야 한다. IBRION=1은 위에서도 말했듯이 어느정도 구조가 minimum 근처에 있을 때 사용하면 좋고, IBRION=3은 그렇지 않을 때 사용하면 좋다. 한편 IBRION=1 또는 3을 사용하게 되면 POTIM을 적절히 설정해 줄 필요가 있는데, POTIM을 작게 할수록 최적화가 보수적으로 일어나기 때문에 느리지만 안정적이게 되고, 크게 할수록 과감해지기 때문에 빠르지만 불안정할 수 있다. 따라서 band relaxation을 처음 시작할 때, relaxed band의 구조에 대한 사전 정보가 전혀 없고 이미지들의 NEB force가 굉장히 높다면 IBRION=3을 사용하고, POTIM은 0.01~0.1 정도로 작은 값을 사용하는 것을 고려해볼 수 있다. 이후 각 이미지의 NEB force가 1 ~ 0.5 eV/A 정도로 줄어든다면 POTIM을 늘리거나 IBRION=1을 사용해볼 수도 있다. 물론 처음에 사용한 보수적인 설정을 계속해서 사용할 수도 있고, 처음부터 과감한 설정을 해서 잘 수렴시킬 수도 있다. 다만 최대한 효율적으로 수렴시키고 싶다면 optimizer를 이와 같이 조절해볼 수 있고, 간혹 optimizer의 설정이 수렴성에 영향을 크게 주는 경우도 있는 것 같다. 그런 의미에서 다양한 optimizer 설정을 테스트해보는 것이 좋다고 생각된다.

 

4.3. Initial trial images

NEB 계산에서 initial image를 잘 설정해주는 게 좋은데, 보통은 원자들의 좌표를 linear interpolation 해서 정해주는 경우가 많다. 이런 initial image를 잘 정하기 위한 알고리즘도 제안된 바가 있는데, Image Dependent Pair Potential (IDPP, https://doi.org/10.1063/1.4878664) 라는 것이 있다. 이건 원자의 좌표가 아니라 원자간 결합 길이를 interpolate하는 것이다. 이 방법은 ASE와 pymatgen-diffusion-analysis 패키지에 구현이 되어 있는데, 나는 ASE를 사용해서 해 본 적은 없고 후자를 사용해 본 적은 있다. 물론 당연히 이게 원자 좌표를 interpolate하는 것보다 무조건 좋으리란 법은 없으니 상황에 맞게 써보면 될 것이다. 다른 알고리즘으로는 Ceder 그룹에서 나온 PathFinder (https://doi.org/10.1063/1.4960790) 가 있는데 이건 ionic material에서 이온의 이동을 모사할 때 사용할 만한 알고리즘으로 제시된 거라 일반적으로도 쓸만한 건지는 모르겠다.

 

Henkelman group에서는 NEB 계산에 특화된 optimizer들을 개발하여 배포하였다. 이에 대한 설명은 [6]에서 확인할 수 있다. Henkelman 교수는 주로 IOPT=3과 IOPT=1을 많이 언급하고 있는 것 같다. 즉 IOPT=3으로 preconverge시킨 후 IOPT=1로 최종 수렴시키는 방식이다. 건드려볼만한 parameter로는 IOPT=3에서는 TIMESTEP, IOPT=1에서는 INVCURV가 있다. TIMESTEP은 POTIM과 비슷한 역할을 하고, 그래서 작을수록 보수적으로 최적화가 된다. INVCURV는 initial curvature인데, 이것의 의미를 이해하려면 사실 IOPT=1에서 사용하는 최적화 알고리즘을 어느정도 이해할 필요가 있다. 결론만 말하자면 이것 역시 작을수록 보수적이게 된다. Henkelman 교수는 보수적인 설정으로 0.001 ~ 0.002 정도의 작은 숫자를 시도해 보라고 언급한 바 있다. 추가로 IOPT=7도 종종 언급되는 것 같은데, 한번 시도해볼 만한 optimizer라고 생각된다.

한편 IBRION=3이나 IOPT=3 or 7은 first order optimizer로서 아주 정확한 force, 즉 EDIFFG가 아주 낮지 않더라도 잘 작동하는 편이지만 IBRION=1이나 IOPT=1은 second order optimizer, 즉 curvature 정보까지 사용하는 optimizer이므로 force가 정확해야 할 필요가 있다. 따라서 EDIFFG를 낮게 설정해주는 것이 좋다.

 

5. 참조한 사이트 및 유용한 사이트 목록 (vaspwiki는 사이트 목록에서 제외하였다)

[1] https://youtube.com/playlist?list=PLjv4kT_-jnZ_NNQ8LyLw96LD4-PM2p6LS 

재료전산모사

[2] https://www.smcm.iqfr.csic.es/docs/vasp/node81.html

Frequently used settings in the INCAR file

[3] http://kitchingroup.cheme.cmu.edu/dft-book/dft.html

Modeling materials using density functional theory

[4] https://theory.cm.utexas.edu/henkelman/forum/

UT theoretical chemistry code forum - Index page

[5] https://sites.tufts.edu/andrewrosen/density-functional-theory/vasp/

VASP | Rosen Review

[6] https://theory.cm.utexas.edu/vtsttools/optimizers.html

FORCE BASED OPTIMIZERS — Transition State Tools for VASP

[7] https://github.com/materialsproject/custodian/blob/master/custodian/vasp/handlers.py

[8] https://github.com/materialsproject/pymatgen/blob/master/pymatgen/io/vasp/sets.py