01. SIESTA 첫 계산 — 1D carbon chain SCF

첫 번째 실습으로 1D carbon chain(cumulene)의 self-consistent field(SCF) 계산을 수행한다. 이 장에서 만드는 4원자 cell은 튜토리얼 전체를 관통하는 예제 시스템이며, 뒤 챕터에서 전극(electrode)과 device의 재료가 된다. fdf 입력 형식, 구조 정의, 기저와 격자, k-점 샘플링, SCF 수렴 기준을 하나씩 짚는다.

학습 목표

fdf 입력 형식의 기본 문법(키워드-값, %block, %include, 단위 표기)을 익힌다
격자·원자 좌표·화학종을 fdf로 정의한다
LCAO 기저(DZP)와 MeshCutoff의 의미를 구분해서 이해한다
1D 시스템의 k-점 샘플링 원칙(transport 축만 조밀하게)을 적용한다
SCF 출력(dDmax, Eharris, E_KS, Fermi energy)을 읽고 수렴을 판정한다

배경

fdf 입력 형식

SIESTA의 입력은 fdf(flexible data format) 형식의 텍스트 파일이다. 규칙은 단순하다.

키워드-값 쌍: 한 줄에 키워드 값 형태로 쓴다. 예: MeshCutoff 300. Ry
키워드 매칭은 관대하다: 대소문자를 구분하지 않고 ., -, _ 문자를 무시한다. MeshCutoff, Mesh.Cutoff, mesh-cutoff는 모두 같은 키워드다.
물리량에는 단위를 붙인다: 300. Ry, 1.29 Ang, 0.01 eV/Ang처럼 값 뒤에 단위를 명시한다.
블록: 여러 줄 데이터는 %block 이름 ... %endblock 이름으로 감싼다. 격자 벡터, 원자 좌표, k-grid가 대표적이다.
파일 분리: %include 파일명으로 다른 fdf 파일을 불러올 수 있다. 구조를 struct.fdf로 분리하고 본 입력에서 include하는 관례를 뒤 챕터(전극/device)에서 사용한다. 이 장에서는 단일 파일로 간다.
주석: # 뒤는 무시된다.
논리값: T/F (또는 true/false).

예제 시스템 — cumulene carbon chain

탄소 원자가 등간격 $d = 1.29$ Å으로 늘어선 1차원 사슬을 cumulene이라고 부른다(모든 C–C 결합이 동등한 이중결합 성격). 이 튜토리얼의 공통 규격은 다음과 같다.

항목	값
원자 간격	$d = 1.29$ Å (등간격)
unit cell	C 4원자, $c = 4d = 5.16$ Å
transport 축	$z$
vacuum	$x = y = 15$ Å

$x$ , $y$ 방향에 15 Å의 진공을 두는 이유는 주기 경계조건 아래에서 이웃 이미지 사슬 간의 인공적인 상호작용을 끊기 위해서다. LCAO 기저는 원자 반경 바깥에서 정확히 0이 되므로, 궤도가 겹치지 않을 만큼의 진공이면 충분하다.

LCAO 기저와 MeshCutoff — 서로 다른 두 개념

SIESTA는 plane-wave가 아니라 원자에 붙은 국소 궤도(LCAO)를 기저로 쓴다. 기저 크기는 PAO.BasisSize로 정하며, 이 튜토리얼은 DZP(double-zeta polarized)를 표준으로 쓴다. DZP는 각 valence 궤도당 radial 함수 2개에 polarization 궤도를 더한 구성으로, 대부분의 production 계산에서 무난한 선택이다.

MeshCutoff는 기저와 무관하게, 전자 밀도와 포텐셜을 표현하는 실공간 grid의 조밀도를 정한다. plane-wave 코드의 ENCUT과 이름이 비슷해 혼동하기 쉬운데 역할이 다르다 — SIESTA에서 기저 품질은 PAO.BasisSize가, grid 적분 정밀도는 MeshCutoff가 담당한다. 이 튜토리얼은 300 Ry를 표준으로 쓴다.

k-점 샘플링 — 1D는 transport 축만

이 시스템은 $z$ 방향으로만 주기적이고 $x$ , $y$ 는 진공이다. 진공 방향은 밴드 분산이 없으므로 k-점 1개(Γ)면 충분하고, $z$ 방향만 조밀하게 샘플링한다. 금속성 사슬은 Fermi 면(1D에서는 Fermi 점) 근처 샘플링에 민감하므로 넉넉하게 $1 \times 1 \times 64$ 를 쓴다.

SCF 수렴 기준 — 왜 $10^{-8}$ 인가

SCF 반복의 수렴은 밀도 행렬(density matrix, DM)의 반복 간 최대 변화 SCF.DM.Tolerance로 판정한다. 이 튜토리얼은 모든 계산에서 1.0d-8을 쓴다. 디폴트 수준의 느슨한 값( $10^{-4}$ )으로도 SCF는 "수렴"하지만, 그 상태의 총에너지는 마지막 자릿수가 흔들려서 구조 간 에너지 차이 비교를 신뢰할 수 없다. 이 튜토리얼에서는 cumulene과 polyyne의 에너지 차이(원자당 meV 수준)를 비교하는 등 미세한 에너지 스케일을 다루므로, 처음부터 타이트한 기준을 습관으로 삼는다.

입력 파일

작업 디렉토리를 만들고 아래 내용을 input.fdf로 저장한다. 전체 파일은 예제 저장소의 code/ch01-siesta-scf에도 있다.

input.fdf
# 01. 1D carbon chain (cumulene) — first SCF
# 4-atom cell, d(C-C) = 1.29 Ang, c = 5.16 Ang, 15 Ang vacuum
# 필요 파일: C.psml (PseudoDojo, PBE, scalar-relativistic)

SystemName   Cumulene C4 chain
SystemLabel  siesta

NumberOfAtoms    4
NumberOfSpecies  1
%block ChemicalSpeciesLabel
  1  6  C
%endblock ChemicalSpeciesLabel

LatticeConstant 1.0 Ang
%block LatticeVectors
  15.000000   0.000000   0.000000
   0.000000  15.000000   0.000000
   0.000000   0.000000   5.160000
%endblock LatticeVectors

AtomicCoordinatesFormat Ang
%block AtomicCoordinatesAndAtomicSpecies
   7.500000   7.500000   0.000000  1
   7.500000   7.500000   1.290000  1
   7.500000   7.500000   2.580000  1
   7.500000   7.500000   3.870000  1
%endblock AtomicCoordinatesAndAtomicSpecies

# --- 기저 / grid ---
PAO.BasisSize  DZP
MeshCutoff     300. Ry

# --- 교환-상관 ---
XC.Functional  GGA
XC.Authors     PBE

# --- k-점: 1D chain, transport 축 z ---
%block kgrid_Monkhorst_Pack
   1   0   0    0.0
   0   1   0    0.0
   0   0  64    0.0
%endblock kgrid_Monkhorst_Pack

# --- SCF ---
MaxSCFIterations   300
SCF.DM.Tolerance   1.0d-8
SCF.Mixer.Weight   0.2
SCF.Mixer.History  15
DM.UseSaveDM       T

# --- 출력 ---
WriteForces        T

라인 해설

시스템 정의

SystemName — 사람이 읽는 설명 문자열. 출력에 그대로 찍히며 계산에는 영향이 없다.
SystemLabel — 모든 출력 파일의 접두어. 여기서는 siesta이므로 siesta.DM, siesta.XV, siesta.EIG 같은 파일이 생성된다.
%block ChemicalSpeciesLabel — 화학종 정의. 각 줄은 종 번호, 원자번호, 라벨 순이다. 라벨 C는 pseudopotential 파일 이름(C.psml)과 연결된다.

구조 정의

LatticeConstant — 격자 벡터의 스케일 인자. 1.0 Ang으로 두면 LatticeVectors의 숫자가 그대로 Å 단위가 된다.
%block LatticeVectors — 3개의 격자 벡터를 행 단위로 쓴다. $x = y = 15$ Å 진공, $z$ 방향 $c = 5.16$ Å.
AtomicCoordinatesFormat Ang — 원자 좌표를 Å 단위 Cartesian으로 쓴다는 선언. Fractional(격자 벡터 기준 분율 좌표)도 자주 쓴다.
%block AtomicCoordinatesAndAtomicSpecies — 각 줄이 x y z 종번호. 4개 원자가 $z = 0, 1.29, 2.58, 3.87$ Å에 등간격으로 놓이고, cell 중앙( $x = y = 7.5$ Å)에 배치된다.

기저·grid·XC·k-점 — 배경 절에서 설명한 값 그대로다. k-grid 블록은 $3 \times 3$ 정수 행렬과 각 행의 shift(마지막 열)로 구성되며, 대각 성분이 각 방향 분할 수다.

SCF 설정

MaxSCFIterations 300 — SCF 최대 반복 수. 수렴 전에 도달하면 계산이 멈추므로 여유 있게 잡는다.
SCF.DM.Tolerance 1.0d-8 — 수렴 기준(배경 절 참조).
SCF.Mixer.Weight 0.2, SCF.Mixer.History 15 — 밀도 혼합(Pulay mixing) 파라미터. 새 밀도를 20%씩 섞고 이전 15 스텝의 이력을 사용한다. 수렴이 진동하면 weight를 줄이는 것이 첫 번째 처방이다.
DM.UseSaveDM T — 디렉토리에 siesta.DM이 있으면 그것을 초기 밀도로 재사용한다. 같은 시스템을 반복 계산할 때 SCF가 크게 빨라진다(restart의 기본기).

SystemLabel이 출력 파일 이름을 결정한다

모든 출력 파일의 접두어는 SystemLabel 값이다. 이 값이 계산마다 뒤죽박죽이면 후처리 스크립트가 파일을 못 찾는 사고가 반복된다. 이 튜토리얼은 관례를 고정한다 — 단독 계산은 siesta, 전극 계산은 Electrode(챕터 05), device 계산은 trans(챕터 06 이후). 뒤 챕터의 TBtrans·sisl 후처리는 전부 이 관례를 전제로 한다.

pseudopotential 배치

설치 장에서 내려받은 C.psml을 작업 디렉토리(input.fdf와 같은 위치)에 복사한다. SIESTA는 ChemicalSpeciesLabel의 라벨과 같은 이름의 pseudopotential 파일을 작업 디렉토리에서 찾는다. 파일이 없으면 즉시 에러로 종료되므로, 실행 전에 다음 두 파일이 있는지 확인한다.

ls
# C.psml  input.fdf

실행

직렬 실행:

siesta < input.fdf > siesta.out

MPI 병렬 빌드라면:

mpirun -np 4 siesta < input.fdf > siesta.out

이 크기(4원자, DZP)의 계산은 노트북에서 수 분 안에 끝난다. 진행 상황은 다른 터미널에서 tail -f siesta.out으로 지켜볼 수 있다.

출력 분석

siesta.out은 크게 헤더(버전·입력 echo) → 기저 생성 로그 → SCF 사이클 → 최종 에너지·힘 순으로 구성된다.

SCF 사이클

파일 중반에서 다음 형태의 표를 찾는다(수치와 열 구성은 버전·빌드에 따라 조금 다르다 — 아래는 예시 발췌).

        iscf     Eharris(eV)        E_KS(eV)     FreeEng(eV)     dDmax    Ef(eV)  dHmax(eV)
   scf:    1     -621.383914     -598.130561     -598.152804  1.240718 -8.309820  4.879274
   scf:    2     -600.560423     -611.238379     -611.260622  0.421559 -5.410181  1.680392
   ...
   scf:   34     -616.283719     -616.283719     -616.305962  0.000000 -4.591428  0.000001

SCF cycle converged after 34 iterations

각 열의 의미:

Eharris — Harris functional 에너지. 수렴 과정에서 E_KS와 다른 경로로 접근하며, 수렴하면 둘이 일치한다. 두 값의 차이가 SCF 수렴의 감각적인 척도다.
E_KS — Kohn–Sham 총에너지.
FreeEng — 전자 smearing에 의한 엔트로피 항까지 포함한 free energy.
dDmax — 밀도 행렬 원소의 반복 간 최대 변화. 이 값이 SCF.DM.Tolerance( $10^{-8}$ ) 아래로 내려가는 것이 수렴 조건이다.
Ef — 해당 반복에서의 Fermi energy (eV).
dHmax — Hamiltonian 원소의 최대 변화. DM 기준과 함께 보조 수렴 지표로 쓰인다.

수렴 판정에서 확인할 것은 두 가지다 — ① SCF cycle converged 메시지가 있는가, ② 마지막 반복의 dDmax가 설정한 tolerance 이하인가. MaxSCFIterations에 걸려 멈춘 계산은 수렴한 것이 아니므로 결과를 쓰면 안 된다.

총에너지와 Fermi energy

파일 후반의 최종 요약에서:

siesta: E_KS(eV) =             -616.2837

siesta: Final energy (eV):
...
siesta:         Total =    -616.283719
siesta:         Fermi =      -4.591428

총에너지(Total)는 pseudopotential 기준의 값이므로 절대값 자체는 의미가 없고, 같은 조건(기저·grid·k-점) 계산 간의 차이만 물리적 의미를 가진다.
Fermi energy는 뒤 챕터에서 밴드·transmission의 에너지 기준점( $E - E_F$ )으로 계속 쓰인다.

생성된 파일

파일	내용
`siesta.DM`	수렴된 밀도 행렬 — 재계산 시 초기 밀도로 재사용(`DM.UseSaveDM`)
`siesta.XV`	최종 구조(좌표+속도)
`siesta.EIG`	샘플링 k-점의 고유값 — 첫 줄이 Fermi energy
`siesta.FA`	원자별 힘 (`WriteForces T`)

등간격 cumulene은 모든 원자가 대칭적으로 동등하므로 siesta.FA의 힘이 0에 가깝게 나오는 것이 정상이다.

연습문제

k-점 수렴 테스트 — $k_z$ 를 16, 32, 64, 128로 바꿔 총에너지를 표로 정리하라. 인접한 두 설정의 에너지 차이가 원자당 1 meV 아래로 내려가는 지점이 어디인지 찾아라. 금속성 1D 시스템에서 k-수렴이 왜 느린지 Fermi 점 샘플링 관점에서 설명해 보라.
MeshCutoff 수렴 테스트 — MeshCutoff를 200, 250, 300, 400 Ry로 바꿔 총에너지 변화를 확인하라. grid가 성겨지면 원자 위치에 따라 에너지가 흔들리는 eggbox 효과가 나타난다 — 300 Ry가 무난한 타협점임을 수치로 확인해 보라.
기저 크기 비교 — PAO.BasisSize를 SZ, DZ, DZP로 바꿔 총에너지와 Fermi energy를 비교하라. 기저가 커질수록 에너지가 낮아지는 이유(variational principle)를 생각해 보라.

다음 장에서는 이 계산에 밴드 경로와 PDOS 블록을 얹어 cumulene이 왜 금속인지 확인한다 — 챕터 02 — 밴드 구조와 DOS.

학습 목표​

배경​

fdf 입력 형식​

예제 시스템 — cumulene carbon chain​

LCAO 기저와 MeshCutoff — 서로 다른 두 개념​

k-점 샘플링 — 1D는 transport 축만​

SCF 수렴 기준 — 왜 10−810^{-8}10−8인가​

입력 파일​

라인 해설​

pseudopotential 배치​

실행​

출력 분석​

SCF 사이클​

총에너지와 Fermi energy​

생성된 파일​

연습문제​