분석 방법

NOTE 32-07-00 · 분석 방법

1. 일반 사항

A. 목적

1. 이 노트는 Cu 표면 벤젠-에탄올 계의 주요 분석 항목과 각 항목이 알려 주는 내용을 다룬다.
2. 분석 항목은 RDF, z 방향 밀도 프로파일, SEI, 흡착 에너지, 계면 장력이다.

B. 적용 범위

1. UROPS 때 후처리에는 integrated_analysis.py (MDAnalysis 기반) 를 썼다. 이 스크립트는 공개하지 않았다.
2. 밀도 프로파일과 첫 층 조성, 블록 평균은 PPPM vs MSM 글의 profiles.py, blocks.py 로 덤프에서 다시 계산할 수 있다.
3. 압력과 계면 장력의 단위는 LAMMPS real 단위를 기준으로 한다. metal 단위 환산은 3.E 에 함께 적는다.

C. 결과 요약

1. SEI 는 절댓값이 있어 음수가 될 수 없다. 방향은 $x^\text{surf} - x^\text{bulk}$ 로 따로 보고한다.
2. UROPS 보고서의 흡착 에너지는 지수형 모형 값이므로 쓰지 않는다. compute group/group 으로 다시 잰다.
3. 이 계는 두 계면이 달라 Irving-Kirkwood 적분 전체를 반으로 나눈 값은 어느 계면의 장력도 아니다.
4. 원자별 응력은 fix ave/chunk 대신 compute reduce/chunk ... sum 으로 빈별 합을 낸다.

2. 준비 정보

A. 필요한 개념

1. 동경 분포 함수(RDF)와 밀도 프로파일의 정의.
2. 표면 과잉량(surface excess)과 계면 장력의 열역학적 정의.
3. compute 와 fix ave/chunk 출력(NOTE 31-07-00).

B. 참조 자료

C. 사용 프로그램

D. 관련 파일

3. 절차

A. 동경 분포 함수 (Radial Distribution Function, RDF)

1. 동경 분포 함수 $g(r)$ 을 정의한다. 기준 입자에서 거리 $r$ 만큼 떨어진 곳에 다른 입자가 있을 확률을 이상 기체(균일 분포)와 비교해 정규화한 양이다.

$$ g_{\alpha\beta}(r) = \frac{1}{\rho_\beta} \left\langle \sum_{i \in \alpha} \sum_{j \in \beta, j \neq i} \frac{\delta(r - r_{ij})}{4\pi r^2} \right\rangle $$

2. 다음 원자 쌍을 분석한다.

3. 벤젠-벤젠 stacking 거리를 ab initio 범위와 비교한다. 벤젠 dimer의 ab initio 계산값은 3.4-3.9 Å (Sherrill 외) 이다.

4. RDF 를 LAMMPS 안에서 바로 계산하려면 다음 명령을 쓴다.

compute rdf_cu_o all rdf 200 12 5 cutoff 12.0     # OPLS-AA: Cu(12)-O(5)
fix     rdf_save all ave/time 100 10 1000 c_rdf_cu_o[*] file rdf_cu_o.dat mode vector

B. 밀도 프로파일 (z 방향)

1. 슬랩 시스템에서 가장 중요한 분석량인 z 방향(표면 수직)의 분자별 밀도 분포를 정의한다.

$$ \rho_\alpha(z) = \frac{\langle N_\alpha(z, z + \Delta z) \rangle}{L_x L_y \Delta z} \cdot m_\alpha $$

2. 여기서 $m_\alpha$ 는 종 $\alpha$ 분자의 질량이다.
3. 예상되는 프로파일 모양과 비교한다.

4. LAMMPS 내부에서 밀도 프로파일을 계산하려면 다음 명령을 쓴다.

compute     z_bins all chunk/atom bin/1d z lower 1.0 units box
fix         density_save all ave/chunk 100 10 1000 z_bins density/mass file density.dat

5. 그룹별 밀도를 따로 보려면 group을 먼저 정의하고 chunk를 그 그룹에 적용한다.

compute     z_bins_benzene benzene chunk/atom bin/1d z lower 1.0 units box
fix         density_benzene benzene ave/chunk 100 10 1000 z_bins_benzene density/mass file density_benzene.dat

C. 분리 효율 지수 (Separation Efficiency Index, SEI)

1. SEI 를 정의한다. SEI는 표면 1차 흡착층에서 두 분자 종이 얼마나 갈라져 있는지를 수치로 나타낸다.

$$ \text{SEI} = \frac{\left| x_{\text{benzene}}^{\text{surf}} - x_{\text{benzene}}^{\text{bulk}} \right|}{1 - x_{\text{benzene}}^{\text{bulk}}} $$

2. 각 기호를 정한다.
   1. $x_{\text{benzene}}^{\text{surf}}$: 표면 1차 흡착층의 벤젠 분율
   2. $x_{\text{benzene}}^{\text{bulk}}$: 벌크 영역의 벤젠 분율
3. 값을 다음 기준으로 해석한다.
   1. SEI = 0: 표면 조성 = 벌크 조성 (분리 없음, 완전 혼합)
   2. SEI = 1: 표면에 벤젠만 흡착 (완전 분리)

4. 방향을 보려면 절댓값 없이 $x^\text{surf} - x^\text{bulk}$ 를 같이 보고한다.
5. Cu 슬랩 맨 위의 z 좌표 $z_{\text{Cu,max}}$ 를 구한다.
6. 표면 영역을 $z_{\text{Cu,max}} < z < z_{\text{Cu,max}} + 5$ Å 로 정한다.
7. 벌크 영역을 $z > z_{\text{Cu,max}} + 10$ Å (또는 박스 위쪽의 일정 거리)로 정한다.
8. 각 영역에서 분자 종별 분자 수(또는 질량)를 평균한다.
9. 분율을 구하고 그 차이를 계산한다.

10. SEI 는 블록 평균과 그 표준오차를 함께 낸다(4.A).

D. 흡착 에너지 (Adsorption Energy)

1. 분자 $\alpha$ 의 흡착 에너지를 다음과 같이 정의한다.

$$ \Delta E_{\text{ads}}^\alpha = E_{\text{surf+mol}} - E_{\text{surf}} - E_{\text{mol}} $$

2. 각 항을 정한다.
   1. $E_{\text{surf+mol}}$: 분자가 표면에 흡착된 시스템의 에너지
   2. $E_{\text{surf}}$: 같은 표면만 있을 때의 에너지
   3. $E_{\text{mol}}$: 같은 분자가 진공(vacuum)에 홀로 있을 때의 에너지
3. production 궤적에서 다음 양으로 흡착 에너지를 근사한다.

$$ E_{\text{ads}}^\alpha \approx \langle E_{\text{Cu-}\alpha}^{\text{LJ+Coul}} \rangle_{\text{surf}} - \langle E_{\text{Cu-}\alpha}^{\text{LJ+Coul}} \rangle_{\text{far}} $$

4. LAMMPS 에서 group-group 상호작용 에너지를 다음 명령으로 계산한다.

group benzene type 1 2
group copper  type 12

# 두 그룹 간 비결합 상호작용 에너지: compute ID 그룹1 group/group 그룹2
compute     cu_benzene copper group/group benzene kspace yes
thermo_style custom step temp c_cu_benzene

5. 표면 가까이 있는 벤젠(1차 흡착층)과 멀리 있는 벤젠(벌크)의 평균값을 비교해 흡착 에너지로 바꾼다.
6. 결과는 분자 하나당으로 나눠 보고한다.

7. UROPS 보고서의 흡착 에너지는 위의 compute group/group 으로 다시 잰다.

E. 계면 장력 (Interfacial Tension): Irving-Kirkwood 방법

1. 계면 장력 $\gamma$ 를 응력 텐서(stress tensor)의 z 방향 성분과 가로 방향 성분의 차이로 정의한다 (Irving & Kirkwood 1950).

$$ \gamma = \frac{1}{2} \int_{-\infty}^{+\infty} \left[ P_{zz}(z) - \frac{P_{xx}(z) + P_{yy}(z)}{2} \right] dz $$

2. 여기서 $P_{xx}, P_{yy}, P_{zz}$ 는 각 위치의 압력 텐서 성분이다.

3. 적분 구간을 계면 하나만 포함하도록 자른다.

4. LAMMPS 에서 다음 명령으로 빈별 응력 합을 계산한다.

compute   stress_atom all stress/atom NULL pair bond angle dihedral improper kspace fix
compute   z_bins all chunk/atom bin/1d z lower 1.0 units box
compute   stress_sum all reduce/chunk z_bins sum c_stress_atom[1] c_stress_atom[2] c_stress_atom[3]
fix       stress_save all ave/time 100 10 1000 c_stress_sum[*] file stress_profile.dat mode vector

5. 빈 $k$ 의 압력 성분은 후처리에서 $P_{\alpha\alpha}(z_k) = -\sum_{i \in k} s_{i,\alpha\alpha} / (L_x L_y \Delta z)$ 로 계산한다.

6. 단위를 변환한다. LAMMPS real 단위에서 압력은 atm, 길이는 Å 이므로 $\gamma$ 는 atm·Å 로 나온다.
7. 1 atm·Å = 101325 Pa × 10⁻¹⁰ m = 1.01325 × 10⁻⁵ N/m 이므로 다음 식으로 mJ/m² 로 바꾼다.

$$ \gamma\,[\text{mJ/m}^2] = \gamma\,[\text{atm·Å}] \times 1.01325 \times 10^{-2} $$

8. metal 단위(압력 bar)라면 1 bar·Å = 10⁻⁵ N/m = 0.01 mJ/m² 로 바꾼다.
9. 이 시스템에서 예상되는 값과 비교한다.

F. 통합 분석 스크립트 실행

1. integrated_analysis.py 하나로 위의 분석을 모두 돌린다.
2. Production 단계만 분석하려면 다음 명령을 쓴다.

# Production 단계 분석
python integrated_analysis.py \
    --topology trappe.data \
    --trajectories 05_production.lammpstrj \
    --stages production \
    --log lammps_run.log \
    --output analysis_trappe_msm

3. heating, equilibration, production 을 모두 분석하려면 다음 명령을 쓴다.

# 다단계 분석 (heating, equilibration, production 모두)
python integrated_analysis.py \
    --topology trappe.data \
    --trajectories 03_heat.lammpstrj 04_eq.lammpstrj 05_production.lammpstrj \
    --stages heating equilibration production \
    --log lammps_run.log \
    --output analysis_trappe_msm

4. 스크립트가 낸 흡착 에너지는 쓰지 않는다.

4. 시험 및 검사

A. 분석 결과의 통계적 신뢰성

1. 자기 상관 시간을 확인한다. $\tau_A$ 가 production 시간보다 충분히 짧아야 한다. (Allen & Tildesley 2017, Frenkel & Smit 2002)
2. 블록 평균(Block average)을 낸다. production 을 5-10 개 블록으로 나눠 블록 사이의 분산을 계산한다.
3. 수렴 그래프를 그린다. 누적 평균이 시간에 따라 수렴하는지 확인한다.