본문으로 건너뛰기

03. VASP으로 구조 최적화

transport 계산의 입력이 되는 구조는 최적화(relaxation)를 거쳐 확정해야 한다. 이 튜토리얼은 구조 확정 단계에 plane-wave 코드인 VASP을 사용하고, 그 결과를 SIESTA로 이관해 교차 검증하는 흐름을 표준으로 삼는다. VASP은 뒤에서 DFT 준위 오차를 보정할 때 hybrid functional reference로도 다시 등장한다(챕터 11 예고).

이 교차 검증 흐름은 carbon chain에 국한되지 않고, 이후 다른 시스템(예: 금속 전극 분자 접합)에서도 그대로 쓰는 표준 절차다. 요약하면:

  1. 격자·좌표 이관 — CONTCAR의 격자 벡터와 좌표를 단위·형식만 바꿔 fdf로 옮긴다(아래 스니펫). 변환 결과의 결합 길이를 원본과 대조한다.
  2. 밴드 개형 대조 — 같은 functional로 두 코드의 밴드 분산·gap 유무를 겹쳐 확인한다.
  3. 상대 에너지 순서 대조 — 구조 A/B 에너지 차이의 부호와 크기를 비교한다. 총에너지 절대값 비교는 금지다 — pseudopotential과 에너지 기준이 코드마다 다르다.

학습 목표​

  • transport 워크플로우에서 VASP의 역할(구조 확정, reference 전자구조)을 이해한다
  • 4대 입력 파일(INCAR/POSCAR/KPOINTS/POTCAR)을 작성하고 POTCAR를 합성한다
  • carbon chain relax를 실행하고 OSZICAR/OUTCAR/CONTCAR를 읽는다
  • 힘 기준(∣F∣<0.01|F| < 0.01 eV/Å) 수렴 판정을 적용한다
  • relaxed 구조를 SIESTA fdf 형식으로 변환하고 두 코드를 교차 검증한다

배경​

왜 구조 최적화에 VASP인가​

SIESTA로도 구조 최적화는 가능하다. 그럼에도 이 튜토리얼이 VASP을 쓰는 이유는 방법론적 독립성 때문이다.

  • plane-wave 기저는 ENCUT 하나로 체계적으로 수렴시킬 수 있고, 기저가 원자에 붙어 있지 않아 LCAO에서 신경 써야 하는 기저 최적화 문제에서 자유롭다.
  • 구조·격자상수를 서로 다른 기저의 두 코드로 계산해 일치를 확인하면, 어느 한쪽 기저의 아티팩트가 아니라는 확신을 얻는다. 같은 PBE functional을 쓰는 한 두 코드의 구조는 거의 일치해야 정상이다.
  • 뒤 챕터에서 hybrid functional 계산(준위 보정의 reference)도 VASP이 담당하므로, 구조 단계부터 같은 코드로 일관성을 유지한다.

4대 입력 파일​

VASP은 고정된 이름의 파일 4개를 작업 디렉토리에서 읽는다.

파일내용
INCAR계산 파라미터 (태그 = 값)
POSCAR격자 벡터와 원자 좌표
KPOINTSk-점 샘플링
POTCARPAW pseudopotential (라이선스 자산)

POTCAR 합성​

POTCAR는 원소별 PAW 파일을 POSCAR의 원소 순서 그대로 이어붙여 만든다. 라이선스 보유 기관이 제공하는 potpaw 데이터베이스 경로를 POTCAR_DIR이라 하면:

cat $POTCAR_DIR/C/POTCAR > POTCAR

이 예제는 탄소 한 종류라 한 줄이면 되지만, 다원소 시스템에서는 POSCAR 6행(원소 목록)의 순서와 cat 순서가 어긋나면 원자에 엉뚱한 pseudopotential이 붙은 채 조용히 계산이 돌아간다. 합성 후 grep TITEL POTCAR로 원소 순서를 반드시 확인하는 습관을 들인다.

POTCAR는 배포 금지

POTCAR는 VASP 라이선스에 묶인 자산이다. 이 사이트의 예제 저장소에는 POTCAR가 포함되어 있지 않으며, 공개 저장소에 POTCAR를 올리는 것은 라이선스 위반이다. 각자 기관의 potpaw 데이터베이스에서 합성해 사용한다.

입력 파일​

예제는 챕터 01·02와 같은 4원자 carbon chain cell이다. 전체 파일은 예제 저장소의 code/ch03-vasp에 있다.

한 가지 중요한 설계가 초기 구조에 들어 있다. 완전 등간격 cumulene은 모든 원자가 대칭적으로 동등해서 힘이 정확히 0이고, relax를 돌려도 첫 스텝에서 그대로 끝날 수 있다. 이는 Peierls 왜곡 방향의 불안정 모드를 수치적으로 탐색하지 못하는 대칭 보존 stationary point다. 그래서 초기 구조에 Peierls 왜곡 방향으로 작은 교란(결합 길이 ±0.02 Å)을 주어, relax가 더 낮은 에너지의 BLA 구조를 탐색하게 한다.

POSCAR
C4 chain (cumulene + 0.02 Ang BLA perturbation)
1.0
   15.0000000000    0.0000000000    0.0000000000
    0.0000000000   15.0000000000    0.0000000000
    0.0000000000    0.0000000000    5.1600000000
C
4
Direct
  0.5000000000  0.5000000000  0.0000000000
  0.5000000000  0.5000000000  0.2461240310
  0.5000000000  0.5000000000  0.5000000000
  0.5000000000  0.5000000000  0.7461240310
  • 1행은 주석, 2행은 전체 스케일 인자.
  • 3–5행이 격자 벡터 — 챕터 01의 fdf와 같은 cell(15×15×5.1615 \times 15 \times 5.16 Å)이다.
  • 6–7행이 원소 목록과 개수. 이 순서가 POTCAR 합성 순서를 결정한다.
  • Direct는 분율 좌표라는 뜻이다. zz 분율 0.2461, 0.7461은 결합 길이가 1.27/1.31 Å로 교대하는 배치다(등간격 1.29 Å에서 ±0.02 Å 교란).
INCAR
SYSTEM = C4 chain relaxation

PREC   = Accurate
ENCUT  = 500
LREAL  = .FALSE.

ISMEAR = 0
SIGMA  = 0.05

EDIFF  = 1E-8
NELM   = 200

IBRION = 2
ISIF   = 2
NSW    = 100
EDIFFG = -0.01

INCAR 해설​

태그값의미
PRECAccurateFFT grid 등 정밀도 프리셋 — 힘 계산의 기본기
ENCUT500plane-wave cutoff (eV). 탄소 PAW potential이 hard해서 넉넉히 잡는다
LREAL.FALSE.작은 cell은 역공간 projection이 더 정확하다
ISMEAR0Gaussian smearing. 금속/반도체 여부가 relax 중 바뀔 수 있는 시스템에 무난한 선택
SIGMA0.05smearing 폭 (eV). 작게 잡아 smearing에 의한 에너지 오염을 줄인다
EDIFF1E-8전자 SCF 수렴 기준 (eV). 힘의 품질은 SCF 수렴에 직결되므로 타이트하게
NELM200SCF 최대 반복 수
IBRION2conjugate-gradient 이온 최적화
ISIF2원자 위치만 relax, cell은 고정
NSW100이온 스텝 최대 수
EDIFFG-0.01음수 = 힘 기준. 모든 원자의 힘이 $
KPOINTS
Gamma-centered 1x1x32 (transport axis z only)
0
Gamma
 1 1 32
 0 0 0

k-점 원칙은 SIESTA 때와 같다 — 진공 방향(xx, yy)은 1, 주기 방향(zz)만 조밀하게. 4원자 cell 기준 32는 primitive cell 환산 128에 해당하는 충분한 샘플링이다.

격자상수까지 최적화하려면 — ISIF 3의 함정

ISIF = 3은 원자와 cell을 함께 relax하지만, 이 시스템에는 그대로 쓰기 어렵다. 응력(stress)이 세 방향 모두에 걸리므로 진공 방향(xx, yy)의 격자까지 변하고, 진공 크기는 물리적으로 최적화할 자유도가 아니다. 특정 축만 골라 relax하는 selective lattice relaxation은 표준 VASP의 기본 기능이 아니다. 따라서 1D chain의 격자상수 cc를 정하는 방법은 두 가지다.

  1. cc 수동 스캔(권장) — cc를 5.00–5.30 Å 범위에서 몇 점 잡고 각각 ISIF = 2 relax를 돌려 E(c)E(c) 곡선의 최소를 찾는다(연습문제 1).
  2. cell 응력을 그대로 두고 ISIF = 2 결과의 zz 방향 응력을 확인해 cc 보정 방향의 참고로만 쓴다.

또한 cell이 변하는 계산에서는 basis set이 cell에 붙어 정의되므로 Pulay stress 오차가 생긴다 — cell 최적화를 하는 경우 ENCUT을 더 높이는 것이 표준 처방이다.

실행​

POTCAR까지 4개 파일을 갖춘 뒤 실행한다. 실행 파일 이름(vasp_std 등)과 병렬 설정은 기관 빌드에 따라 다르다.

mpirun -np 4 vasp_std > vasp.out

표준 출력을 vasp.out으로 리다이렉트해 실행 로그를 남겨 둔다 — VASP은 SIESTA와 달리 로그를 파일로 자동 저장하지 않으므로, 리다이렉트하지 않으면 실행 기록이 사라진다.

출력 분석​

OSZICAR — 이온 스텝별 에너지​

       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.315293846529E+02   -0.31529E+02   ...
...
   1 F= -.31581462E+02 E0= -.31581523E+02  d E =-.315815E+02
DAV:   1    ...
...
   2 F= -.31581891E+02 E0= -.31581952E+02  d E =-.428999E-02
   ...
  • DAV: 줄은 전자 SCF 반복, 맨 앞에 스텝 번호가 붙은 F= 줄이 이온 스텝 하나의 완료다.
  • F는 free energy, E0은 smearing을 0으로 외삽한 에너지다. ISMEAR = 0을 쓸 때 구조 간 에너지 비교는 E0으로 한다.
  • d E가 이온 스텝마다 줄어들며 0으로 접근하는지 본다.

OUTCAR — 힘과 수렴 메시지​

각 이온 스텝의 원자별 힘은 OUTCAR의 TOTAL-FORCE 표에 있다.

grep -A 8 "TOTAL-FORCE" OUTCAR | tail -9
 POSITION                                       TOTAL-FORCE (eV/Angst)
 -----------------------------------------------------------------------------------
      7.50000      7.50000      0.00000       0.000000    0.000000    0.004213
      ...

수렴 판정은 두 가지로 한다 — ① 마지막 스텝의 모든 힘 성분이 0.01 eV/Å보다 작은가, ② vasp.out(또는 OUTCAR)에 reached required accuracy 메시지가 있는가. NSW를 다 쓰고 멈춘 계산은 수렴이 아니다.

CONTCAR — 최종 구조​

CONTCAR는 POSCAR와 같은 형식으로 기록된 최종(또는 마지막 스텝) 구조다. 이 예제에서 물리적으로 확인할 것은 최종 BLA다 — 분율 좌표에서 결합 길이를 계산해 보면, PBE는 약한 결합 교대를 남긴다(PBE가 BLA를 과소평가하는 것은 알려진 경향이고, hybrid functional은 더 큰 BLA를 준다 — 챕터 11에서 다시 만난다).

SIESTA로 구조 이관​

relaxed 구조를 transport 워크플로우(SIESTA)로 가져간다. 단위·형식 변환이 핵심이다 — CONTCAR는 스케일 인자가 곱해지는 분율(Direct) 좌표, 우리의 fdf 관례는 Å 단위 Cartesian 좌표다.

import numpy as np

with open("CONTCAR") as f:
    lines = f.readlines()

scale = float(lines[1].split()[0])
cell = scale * np.array(
    [[float(x) for x in lines[i].split()[:3]] for i in (2, 3, 4)]
)
symbols = lines[5].split()
counts = [int(x) for x in lines[6].split()]
natoms = sum(counts)

i = 7
if lines[i].strip().lower().startswith("s"):      # Selective dynamics 줄 건너뜀
    i += 1
cartesian = lines[i].strip().lower()[0] in "ck"   # Cartesian/Direct 판별
pos = np.array(
    [[float(x) for x in lines[i + 1 + n].split()[:3]] for n in range(natoms)]
)
if not cartesian:
    pos = pos @ cell                              # 분율 -> Angstrom Cartesian

species = []
for s, c in zip(symbols, counts):
    species += [s] * c
spmap = {s: j + 1 for j, s in enumerate(dict.fromkeys(species))}

print("LatticeConstant 1.0 Ang")
print("%block LatticeVectors")
for v in cell:
    print(f"  {v[0]:14.8f} {v[1]:14.8f} {v[2]:14.8f}")
print("%endblock LatticeVectors")
print()
print("AtomicCoordinatesFormat Ang")
print("%block AtomicCoordinatesAndAtomicSpecies")
for r, s in zip(pos, species):
    print(f"  {r[0]:14.8f} {r[1]:14.8f} {r[2]:14.8f}  {spmap[s]}")
print("%endblock AtomicCoordinatesAndAtomicSpecies")

출력된 블록을 fdf의 해당 블록과 교체하면 된다. sisl로도 같은 변환이 가능하다(sisl.get_sile("CONTCAR").read_geometry()로 읽어 fdf로 저장). 어느 경로를 쓰든 변환 결과의 결합 길이를 원본 CONTCAR와 대조하는 검증을 생략하지 않는다.

교차 검증​

이관한 구조로 SIESTA SCF(챕터 01)와 밴드(챕터 02)를 다시 돌려 두 코드를 대조한다.

  • 비교해도 되는 것 — 평형 격자상수·결합 길이, 밴드 개형(분산, gap 유무), 상대 에너지(구조 A와 B의 에너지 차이). 같은 PBE라면 기저 차이(LCAO vs plane-wave)에도 불구하고 서로 근접해야 하며, 크게 어긋나면 어느 한쪽의 수렴 파라미터(기저, ENCUT, k-점)를 의심한다.
  • 비교하면 안 되는 것 — 총에너지 절대값. pseudopotential과 에너지 기준이 코드마다 달라 절대값 비교는 무의미하다.

연습문제​

  1. 격자상수 cc 스캔 — 등간격 구조로 cc를 5.00, 5.08, 5.16, 5.24, 5.32 Å로 바꿔(각 POSCAR에서 3–5행과 분율 좌표는 그대로, cc만 수정) single-point 에너지 E0을 구하고 E(c)E(c)에 포물선을 맞춰 평형 cc를 구하라. 이 튜토리얼의 공통값 5.16 Å과 비교해 보라.
  2. 초기 교란 의존성 — 초기 BLA 교란을 0.00, 0.02, 0.05 Å로 바꿔 relax를 돌리고 최종 BLA를 비교하라. 교란 0.00에서 relax가 왜 등간격에 그대로 머무는지(대칭 논리), 그리고 교란을 준 경우의 최종 BLA가 교란 크기와 무관하게 같은 값으로 수렴하는지 확인하라.
  3. 교차 검증 실습 — relax된 CONTCAR를 위 스니펫으로 fdf에 이관해 SIESTA 밴드를 계산하고, 챕터 02의 등간격 cumulene 밴드와 겹쳐 그려라. BLA가 생기면서 Γ 부근 Fermi 교차점에 gap이 열렸는지, gap 크기가 얼마인지 확인하라.

구조가 확정되었으니 transport로 넘어갈 준비가 됐다. 다음 장부터 NEGF 형식론과 전극 계산이 시작된다.