본문으로 건너뛰기

IP/EA 계산 — ΔSCF와 level alignment의 재료

챕터 11의 DFT+Σ 보정에서 gas-phase 성분 Σ1\Sigma^1의 크기는 고립 분자의 IP(ionization potential)/EA(electron affinity)와 DFT eigenvalue의 차이에서 온다. 그 IP/EA를 실제로 어떻게 계산하는지가 이 장의 주제다. 방법 자체는 total-energy 차이(ΔSCF)라는 단순한 정의지만, 주기 경계 코드(VASP)에서 하전 셀을 다루는 실무에는 알려진 함정이 여럿 있다. 정의 → Koopmans 근사와의 차이 → VASP ΔSCF 절차 → 하전 셀 보정의 한계 → 접합 안에서의 image-charge 재규격화 순으로 정리한다.

학습 목표

  • IP/EA를 total-energy 차이로 정의하고, Koopmans 근사(εHOMO-\varepsilon_{\mathrm{HOMO}})와의 차이를 설명할 수 있다.
  • VASP에서 NELECT로 전자 수를 바꿔 ΔSCF IP/EA를 계산하는 절차를 수행할 수 있다.
  • 주기 경계 하전 셀의 인위적 정전 상호작용과 그 보정(box 크기 외삽)의 한계를 안다.
  • gas-phase IP/EA가 접합 안에서 전극 스크리닝으로 재규격화되는 방향과 크기 어림을 설명할 수 있다.

정의 — total-energy 차이가 정의다

NN 전자 계의 IP와 EA는 전자 수가 다른 계들의 바닥상태 총에너지 차이로 정의된다.

IP=E(N1)E(N),EA=E(N)E(N+1)\mathrm{IP} = E(N-1) - E(N), \qquad \mathrm{EA} = E(N) - E(N+1)

두 개(또는 세 개)의 SCF 계산의 총에너지를 빼서 구하므로 ΔSCF라 부른다. 구조를 어느 시점에 고정하느냐에 따라 두 종류가 있다.

  • vertical — 중성 평형 구조를 고정한 채 전자 수만 바꾼다. 광전자 분광의 수직 전이, 그리고 챕터 11의 level alignment 보정에 들어가는 값이 이것이다(접합 안 분자의 구조는 하전 순간에 이완할 시간이 없다).
  • adiabatic — 하전 상태에서 구조를 다시 relax한다. 열역학적 이온화 에너지에 해당한다.

이 장의 목적(준위 보정의 재료)에는 vertical을 쓴다. 어느 쪽을 썼는지 결과 보고에 항상 명시한다.

두 값의 차이가 분자의 fundamental gap이다.

Egapfund=IPEAE_{\mathrm{gap}}^{\mathrm{fund}} = \mathrm{IP} - \mathrm{EA}

챕터 11에서 "gas-phase gap"이라 부른 것이 바로 이 양이고, PBE Kohn–Sham gap이 이보다 크게 작다는 것이 DFT+Σ의 첫 번째 동기였다.

Koopmans 근사와 그 실패

DFT 한 번의 계산에서 IP를 읽는 지름길로 εHOMO-\varepsilon_{\mathrm{HOMO}}(Koopmans 근사)가 자주 쓰인다. 정확한 Kohn–Sham 이론에서는 εHOMO=IP-\varepsilon_{\mathrm{HOMO}} = \mathrm{IP}가 정리로 성립하지만, semi-local functional의 self-interaction 오차는 점유 준위를 위로 밀어 올려 이 관계를 크게 깨뜨린다. benzene이 잘 알려진 예다.

방법benzene IP 추정 (eV, 대략값)비고
실험 (vertical)9.2광전자 분광
ΔSCF (PBE)9.2 부근총에너지 차이 — 오차 상쇄로 실험에 근접
εHOMO-\varepsilon_{\mathrm{HOMO}} (PBE)6 대 초반약 3 eV 과소평가 — self-interaction 오차
εHOMO-\varepsilon_{\mathrm{HOMO}} (hybrid, B3LYP/PBE0)7 대정확한 교환 혼합만큼 개선되지만 여전히 과소

요점 두 가지. 첫째, PBE의 εHOMO-\varepsilon_{\mathrm{HOMO}}는 IP를 수 eV 단위로 과소평가하므로 준위 보정의 기준으로 그대로 쓸 수 없다. 둘째, 같은 PBE라도 ΔSCF는 총에너지 차이라 오차가 상당 부분 상쇄되어 작은 분자에서 실험에 훨씬 가깝다 — DFT+Σ의 gas-phase 성분이 eigenvalue가 아니라 ΔSCF(또는 hybrid eigenvalue를 매개로 한 근사, 챕터 11)를 기준으로 삼는 이유다.

왜 ΔSCF가 되는지도 짚어 둔다. semi-local functional의 총에너지는 정수 전자 수에서 상대적으로 정확하고, eigenvalue가 짊어지는 두 오차 — self-interaction과 derivative discontinuity 누락(챕터 11) — 는 E(N1)E(N-1)E(N)E(N)차이에서 상당 부분 상쇄된다. 즉 ΔSCF는 eigenvalue를 우회해 총에너지라는 functional의 더 신뢰할 수 있는 면만 쓰는 전략이다. 다만 이 상쇄는 보장이 아니라 경험적 경향이므로, 가능한 계에서는 실험 vertical IP와의 대조를 함께 제시한다.

실무 절차 — VASP ΔSCF

고립 분자를 주기 코드에서 다루므로, 주기 이미지 간 상호작용을 억제하는 설정이 우선이다.

  • 큰 box + Γ\Gamma — 분자 가장자리 기준 진공 12 Å 이상을 확보하고 k-점은 Γ\Gamma 하나만 쓴다(고립 분자는 분산이 없어야 정상이다).
  • 전자 수 제어NELECT로 전자 수를 ±1\pm1 바꾼다. 중성 계산의 기본 전자 수는 OUTCAR의 NELECT 값으로 확인하고, 거기서 1을 더하거나 뺀다.
  • spin-polarized 필수 — 닫힌 껍질 분자에서 전자를 하나 빼거나 더하면 홀수 전자 계가 되므로 하전 계산은 ISPIN = 2로 돌리고, 최종 자기모멘트가 기대값(1 μB)인지 확인한다.
  • smearing — 분자 준위는 이산적이므로 ISMEAR = 0에 작은 SIGMA(0.01 eV 수준)를 쓴다.
  • NELECT 외 모든 설정 동일 — 세 계산은 ENCUT·box·SIGMA·pseudopotential이 완전히 같아야 한다. 하나라도 다르면 그 차이가 상쇄되지 않고 IP/EA에 그대로 남는다(연습문제 4).
INCAR (양이온, N-1 전자 — vertical IP용)
SYSTEM = molecule cation, fixed neutral geometry
PREC = Accurate
ENCUT = 500
ISMEAR = 0
SIGMA = 0.01
EDIFF = 1E-8
ISPIN = 2
NELECT = 41 # 중성 42 전자의 예 — 중성 OUTCAR의 NELECT에서 1 감소

절차는 세 번의 single-point 계산이다(구조는 중성 relax 결과로 고정).

단계계산전자 수산출
1중성NN (기본값)E(N)E(N)
2양이온NELECT =N1= N-1E(N1)E(N-1)IP=E(N1)E(N)\mathrm{IP} = E(N-1) - E(N)
3음이온NELECT =N+1= N+1E(N+1)E(N+1)EA=E(N)E(N+1)\mathrm{EA} = E(N) - E(N+1)

세 계산의 POSCAR는 동일 파일(중성 relax의 CONTCAR)을 그대로 복사한다 — vertical 정의가 요구하는 조건이다. 전자 수와 에너지는 출력에서 직접 확인한다.

grep NELECT OUTCAR # 이 계산이 실제로 쓴 전자 수 — 의도값과 대조
grep "E0=" OSZICAR | tail -1 # smearing 외삽 총에너지 E0
grep "mag=" OSZICAR | tail -1 # 최종 자기모멘트 — 하전 계산에서 1 μB 근방인지

에너지 비교는 챕터 03과 같은 이유로 E0(smearing 외삽값)로 한다. 하전 계산이 자기모멘트 0으로 수렴했다면 홀수 전자 계가 잘못된 스핀 상태에 앉은 것이므로, 초기 모멘트를 걸어 doublet 해를 다시 찾는다.

참고 — εHOMO-\varepsilon_{\mathrm{HOMO}}를 비교할 때는 진공 준위 정렬이 선행된다

ΔSCF는 같은 box 안의 총에너지 차이라 기준 문제가 자동으로 상쇄되지만, eigenvalue εHOMO\varepsilon_{\mathrm{HOMO}} 자체를 IP와 비교하려면 기준 정렬이 필요하다. 주기 코드의 eigenvalue는 셀 평균 퍼텐셜 기준이므로, LVTOT로 저장한 LOCPOT의 진공 영역 평균값을 진공 준위로 취해 εHOMOVvac\varepsilon_{\mathrm{HOMO}} - V_{\mathrm{vac}}로 읽어야 실험 IP·다른 코드와 비교 가능한 수치가 된다. 본문 표의 εHOMO-\varepsilon_{\mathrm{HOMO}} 값들은 모두 이 정렬을 거친 값 기준이다.

하전 셀의 함정 — 정직하게 쓰기

주기 경계에서 하전 셀은 그 자체로 잘 정의되지 않는다. 코드는 중성화를 위해 균일 배경 전하(jellium)를 자동으로 깔고, 그 결과 총에너지에는 하전 분자·주기 이미지·배경 전하 사이의 인위적 정전 상호작용이 섞인다. 선도 항은 box 크기 LL에 대해 q2/Lq^2/L 스케일로 감쇠한다(Makov–Payne).

  • box 크기 시리즈가 기본 처방이다. box 한 변을 단계적으로 키워(예: 12 → 16 → 20 Å) EE1/L1/L에 대해 외삽한다. 중성 계산은 box에 거의 무관하지만 하전 계산은 그렇지 않으므로, IP/EA의 box 수렴을 수치로 제시해야 한다.
  • 자동 보정 태그는 확인 후 사용. VASP에는 하전·쌍극자 셀에 대한 보정 기능이 있으나 적용 조건이 계산 종류에 따라 다르므로, 사용한 보정과 그 전제를 매뉴얼에서 확인해 명시한다. 어느 쪽이든 box 시리즈로 수렴을 직접 확인하는 것이 가장 논쟁의 여지가 없는 절차다.
  • 음이온은 더 어렵다. gas-phase EA가 음수인 분자(benzene이 그렇다)는 여분 전자가 분자에 속박되지 않는다. 주기 셀에서는 그 전자가 진공 영역의 인위적 상태로 빠져 box를 키울수록 결과가 표류하므로, EA의 ΔSCF는 속박 음이온이 확인되는 계에서만 신뢰할 수 있다. 이 한계는 정직하게 보고한다.
하전 셀 총에너지의 절대값을 믿지 말 것

하전 셀 계산의 총에너지 절대값은 배경 전하와 평균 퍼텐셜 규약에 의존한다. 의미가 있는 것은 보정·외삽을 거친 에너지 차이뿐이며, box 크기와 보정 방법을 명시하지 않은 IP/EA 수치는 재현할 수 없다.

절차 전체를 한 표로 요약한다.

항목설정/판정근거
box진공 12 Å 이상 + 크기 시리즈 외삽주기 이미지·배경 전하 상호작용 q2/L\sim q^2/L
k-점Γ\Gamma 한 점고립 분자 — 분산이 없어야 정상
구조세 계산 모두 중성 relax 구조 고정vertical 정의
스핀하전 계산 ISPIN = 2, 최종 모멘트 확인홀수 전자 계
smearingISMEAR = 0, SIGMA 0.01 eV 수준, 비교는 E0이산 준위 오염 방지
EA속박 음이온 확인 후에만 신뢰비속박 전자의 box 상태 유입

접합 안에서는 gas-phase 값이 그대로 쓰이지 않는다

image charge 준위 재규격화

그림 1. gas-phase HOMO/LUMO가 접합 안에서 전극 스크리닝으로 재규격화되는 그림 — 점유 준위는 올라가고 비점유 준위는 내려가 gap이 2Σimage2\Sigma_{image}만큼 줄어든다. (개념도)

이렇게 얻은 IP/EA는 진공 중 고립 분자의 값이다. 분자를 금속 전극 사이에 넣으면 전자를 빼고 넣을 때 전극에 유도되는 image charge가 하전 상태를 안정화하므로, 점유 준위는 올라가고 비점유 준위는 내려가 준위들이 gap 중앙 쪽으로 이동한다 — gap이 gas-phase보다 좁아진다. 대표 수치로 graphite 위 benzene의 gap이 10.5 eV에서 7.2 eV로 줄어드는 GW 계산이 있다(Neaton, Hybertsen, Louie 2006). 분자 자체의 궤도 모양은 거의 변하지 않는데 준위 위치만 수 eV 이동하는 — 즉 semi-local functional이 원리적으로 담지 못하는 — 비국소 분극 효과라는 것이 이 논문의 요지다.

크기 어림은 고전 정전기로 충분하다. 평판 금속 전극에서 거리 dd에 있는 점전하의 image-charge 안정화 에너지는

Σ=e216πϵ0d\Sigma = \frac{e^2}{16\pi\epsilon_0 d}

이고, e2/4πϵ0=14.4 eVA˚e^2/4\pi\epsilon_0 = 14.4\ \mathrm{eV\cdot\text{Å}}를 쓰면 Σ3.6/d\Sigma \approx 3.6/d eV(dd는 Å)다. d=3d = 3 Å이면 준위당 약 1.2 eV — 전극이 둘이면 image 급수가 더해져 효과는 더 커진다. 주의점:

  • 이 공식은 평면 금속 표면과 잘 정의된 image plane을 전제한다. image plane 위치의 수백 mÅ 불확실성이 보정값 수백 meV를 좌우하고, 평판이 아닌 전극 기하(1D chain 전극이 극단적 예)에서는 공식 자체가 성립하지 않는다 — 챕터 11의 한계 절에서 본 그대로다.
  • 챕터 11의 표기로 연결하면, 점유 준위 보정 Wocc=+ΣW_{\mathrm{occ}} = +\Sigma(위로), 비점유 준위 보정 Wvirt=ΣW_{\mathrm{virt}} = -\Sigma(아래로)다. gas-phase 성분 Σ1\Sigma^1은 gap을 넓히고 image-charge 성분 Σ2\Sigma^2는 gap을 다시 좁힌다 — 두 성분은 부호가 반대다.
  • 따라서 준위 보정의 전체 구조는 (i) 이 장의 gas-phase ΔSCF가 Σ1\Sigma^1을, (ii) 전극 기하에 맞는 image-charge 모델이 Σ2\Sigma^2를 담당하는 2단 구성이다. 어느 성분에 어떤 근사(box 크기·보정, image plane 위치)를 썼는지 분리해 보고해야 결과를 검증할 수 있다.
파이프라인에서의 위치

이 장의 산출물이 들어가는 자리를 다시 확인해 둔다. 챕터 11의 보정식에서 gas-phase ΔSCF의 IP/EA는 Σ1=(IPεHOMODFT)\Sigma^1 = (-\mathrm{IP} - \varepsilon_{\mathrm{HOMO}}^{\mathrm{DFT}}) 항(비점유 쪽은 EA로 대응)을 결정하고, 이 장 마지막 절의 image-charge 어림이 Σ2\Sigma^2를 결정한다. 즉 IP/EA 계산의 오차 막대가 그대로 보정된 공명 위치의 오차 막대가 되고, gap 안 T(EF)T(E_F)는 그 위치에 자릿수 단위로 민감하다 — box 수렴과 보정 방법 명시가 형식 요건이 아니라 결과의 신뢰 구간 그 자체인 이유다.

연습문제

  1. 작은 닫힌 껍질 분자 하나(예: benzene)로 vertical IP를 ΔSCF로 계산하고, 같은 계산의 εHOMO-\varepsilon_{\mathrm{HOMO}}(PBE)와 비교하라. 두 값의 차이가 본문 표의 경향과 맞는지 확인하라.
  2. 양이온 계산을 box 한 변 12, 16, 20 Å에서 반복해 E(N1)E(N-1)1/L1/L에 대해 그리고, 외삽값과 12 Å 값의 차이(= box 유한 크기 오차)를 eV 단위로 보고하라.
  3. 전극-분자 거리 dd를 2.5, 3.0, 4.0 Å로 가정하고 평판 image-charge 공식으로 준위당 Σ2\Sigma^2를 계산하라. gas-phase 보정 Σ1\Sigma^1이 gap을 넓히는 방향, Σ2\Sigma^2가 좁히는 방향임을 부호와 함께 정리하라.
  4. 세 계산(중성/양이온/음이온)에서 NELECT 외의 어떤 설정이라도 달라지면 ΔSCF가 어떻게 오염되는지 논하라 — 특히 ENCUT, box 크기, SIGMA가 계산마다 다를 때 어떤 항이 상쇄되지 않고 남는지 구체적으로 짚어라.

참고문헌

  • J. B. Neaton, M. S. Hybertsen, and S. G. Louie, "Renormalization of Molecular Electronic Levels at Metal-Molecule Interfaces", Phys. Rev. Lett. 97, 216405 (2006). DOI 10.1103/PhysRevLett.97.216405
  • S. Y. Quek, L. Venkataraman, H. J. Choi, S. G. Louie, M. S. Hybertsen, and J. B. Neaton, "Amine-Gold Linked Single-Molecule Circuits: Experiment and Theory", Nano Lett. 7, 3477 (2007). DOI 10.1021/nl072058i
  • G. Makov and M. C. Payne, "Periodic boundary conditions in ab initio calculations", Phys. Rev. B 51, 4014 (1995). DOI 10.1103/PhysRevB.51.4014