Defect Calculation

日志示例 — 原胞来源：

[2026-03-08 10:00:00] STEP: Primitive Cell Acquisition
  Description: 获取 CdTe 原胞结构
  Source: Materials Project, mp-406
  Space group: F-43m (216)
  Lattice: a=b=c=6.629 Å, α=β=γ=60°
  Atoms: 2 (Cd: 1, Te: 1)
  File: primitive_POSCAR

日志示例 — 电荷态分析：

[2026-03-08 14:00:00] STEP: Charge State Analysis - V_Cd
  Description: 分析 V_Cd 中性缺陷的 EIGENVAL，确定可能的电荷态
  Bulk reference:
    VBM: Band 286-288, E=1.790 eV (3-fold degenerate)
    CBM: Band 289, E=2.358 eV
    Band gap: 0.57 eV
  Defect levels identified:
    Band 281-283: E=1.751 eV, occ_up=0.667, occ_down=0.667
    Criteria: partial occupation + near VBM
    Electrons in defect levels: 3×0.667×2 = 4.0 e⁻
    Empty states in defect levels: 2.0
  Charge states determined: q=0, q=-1, q=-2
  NELECT values: q0=564, q-1=565, q-2=566
  Reasoning: 3-fold defect level near VBM, 4 electrons present,
             can accept 2 more → q=-1, q=-2

Calculation Level Framework

Three accuracy levels are supported. The level determines which functional is used for relaxation and static calculations:

Key differences between PBE and HSE INCAR parameters:

Important notes for HSE calculations:

• IALGO = 38 is incompatible with HSE — must use ALGO = Damped (relax) or ALGO = All (static)
• LREAL = .FALSE. is required for HF accuracy (no real-space projection)
• HSE relaxation is ~10-50× more expensive than PBE
• HSE static is ~10-20× more expensive than PBE
• For LEVEL2, band_gap in formation_config.yaml should use the HSE value (not PBE)
• Chemical potential references (mu_O_ref, etc.) should also be recalculated at the HSE level for LEVEL2/3

AEXX Determination (LEVEL2/3 必须步骤)

When user chooses LEVEL2 or LEVEL3, you MUST ask them one of the following three options:

1. Standard HSE06 — Use default AEXX = 0.25. Simplest, suitable when gap accuracy is not critical.
2. User-specified AEXX — User provides a custom AEXX value (e.g., from literature). Use directly.
3. Fit to experimental band gap — You search literature for the experimental band gap, then fit AEXX by running HSE calculations at multiple AEXX values on the primitive cell.

If user chooses Option 3 (fit to experiment), follow this procedure:

Step A: Find Experimental Band Gap

Use web search tools (Tavily, WebSearch) to find the experimental band gap of the bulk material. Search queries like:

• "{material_name} experimental band gap"
• "{material_name} optical band gap measurement"

Record the source (paper DOI, review article, or textbook) and the gap value. Log to defect_workflow.log:

[YYYY-MM-DD HH:MM:SS] STEP: AEXX Fitting - Experimental Band Gap
  Description: 搜索 {material} 实验带隙
  Source: {DOI or reference}
  Experimental band gap: {value} eV
  Type: {direct/indirect, optical/fundamental}

Step B: Prepare AEXX Fitting Calculations

All fitting calculations go in Bulk/AEXX/. Use the primitive cell (NOT the supercell) with a dense k-mesh for accurate band gap determination.

Directory structure:

Bulk/AEXX/
├── primitive_POSCAR          # Primitive cell (from original source or relaxed)
├── aexx_0.20/                # HSE static at AEXX=0.20
│   ├── POSCAR INCAR KPOINTS POTCAR submit_job
├── aexx_0.25/                # HSE static at AEXX=0.25
│   ├── POSCAR INCAR KPOINTS POTCAR submit_job
├── aexx_0.30/                # HSE static at AEXX=0.30
│   ├── POSCAR INCAR KPOINTS POTCAR submit_job
├── aexx_0.35/                # (optional, if needed for wide-gap materials)
│   ├── ...
├── fit_aexx.py               # Fitting script
└── aexx_fit.log              # Results

INCAR for AEXX fitting (static on primitive cell):

SYSTEM = AEXX fitting - primitive cell

# Electronic
ENCUT  = <same as supercell>
EDIFF  = 1E-6          # Tighter for accurate gap
PREC   = Accurate
NELM   = 200
LREAL  = .FALSE.
ISMEAR = 0
SIGMA  = 0.05
GGA    = PE

# HSE
LHFCALC  = .TRUE.
HFSCREEN = 0.2
AEXX     = <varies: 0.20, 0.25, 0.30, ...>
ALGO     = All
TIME     = 0.4
PRECFOCK = Fast

# Spin
ISPIN  = 2

# Output
LORBIT = 10
LWAVE  = .FALSE.
LCHARG = .FALSE.

# Parallelization
NPAR   = 8

KPOINTS for primitive cell: Use the converged k-mesh from Bulk/ktest/ (see "K-mesh Determination" section). The k-mesh convergence test on the primitive cell is mandatory before running AEXX fitting calculations. Example for Al₂O₃ rhombohedral primitive (10 atoms, after convergence test):

Gamma
0
G
6 6 3
0 0 0

Step C: Extract Band Gaps and Fit

After all AEXX calculations converge, extract band gaps from EIGENVAL files and perform linear regression.

fit_aexx.py:

#!/usr/bin/env python3
"""
Fit AEXX to reproduce experimental band gap.
Reads EIGENVAL from Bulk/AEXX/aexx_*/
Usage: python fit_aexx.py --exp_gap <value> --aexx_dir Bulk/AEXX/
"""
import numpy as np
import os
import sys
import argparse
from datetime import datetime

def get_band_gap_from_eigenval(eigenval_path):
    """Extract band gap from EIGENVAL (VBM to CBM)"""
    with open(eigenval_path) as f:
        lines = f.readlines()
    header = lines[5].split()
    nelect = int(float(header[0]))
    nkpts = int(header[1])
    nbands = int(header[2])

    vbm = -999.0
    cbm = 999.0
    idx = 7
    for k in range(nkpts):
        idx += 1  # blank line
        idx += 1  # k-point header
        for b in range(nbands):
            parts = lines[idx].split()
            e_up = float(parts[1])
            occ_up = float(parts[2])
            if len(parts) > 3:
                e_dn = float(parts[3])
                occ_dn = float(parts[4])
                if occ_up > 0.5:
                    vbm = max(vbm, e_up)
                else:
                    cbm = min(cbm, e_up)
                if occ_dn > 0.5:
                    vbm = max(vbm, e_dn)
                else:
                    cbm = min(cbm, e_dn)
            else:
                if occ_up > 0.5:
                    vbm = max(vbm, e_up)
                else:
                    cbm = min(cbm, e_up)
            idx += 1
    return cbm - vbm, vbm, cbm

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--exp_gap', type=float, required=True, help='Experimental band gap (eV)')
    parser.add_argument('--aexx_dir', default='Bulk/AEXX/', help='Directory with aexx_* subdirs')
    args = parser.parse_args()

    # Collect data points
    aexx_values = []
    gaps = []
    for d in sorted(os.listdir(args.aexx_dir)):
        if d.startswith('aexx_'):
            aexx = float(d.replace('aexx_', ''))
            eigenval = os.path.join(args.aexx_dir, d, 'EIGENVAL')
            if os.path.exists(eigenval):
                gap, vbm, cbm = get_band_gap_from_eigenval(eigenval)
                aexx_values.append(aexx)
                gaps.append(gap)
                print(f"  AEXX={aexx:.2f}: gap={gap:.4f} eV (VBM={vbm:.4f}, CBM={cbm:.4f})")

    if len(aexx_values) < 2:
        print("ERROR: Need at least 2 AEXX data points for fitting")
        sys.exit(1)

    # Linear fit: gap = a * AEXX + b
    aexx_arr = np.array(aexx_values)
    gap_arr = np.array(gaps)
    coeffs = np.polyfit(aexx_arr, gap_arr, 1)
    a, b = coeffs

    # Solve for target gap: AEXX_fit = (exp_gap - b) / a
    aexx_fit = (args.exp_gap - b) / a

    result = f"""
=== AEXX Fitting Result ===
Date: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}
Experimental band gap: {args.exp_gap:.4f} eV

Data points:
{''.join(f'  AEXX={av:.2f}: gap={gv:.4f} eV' + chr(10) for av, gv in zip(aexx_values, gaps))}
Linear fit: gap = {a:.4f} * AEXX + {b:.4f}
R² = {1 - np.sum((gap_arr - np.polyval(coeffs, aexx_arr))**2) / np.sum((gap_arr - np.mean(gap_arr))**2):.6f}

>>> Fitted AEXX = {aexx_fit:.4f} <<<
(to reproduce experimental gap of {args.exp_gap:.4f} eV)

Recommendation:
  Round to AEXX = {round(aexx_fit, 2):.2f}
  Verify: expected gap at AEXX={round(aexx_fit, 2):.2f} = {np.polyval(coeffs, round(aexx_fit, 2)):.4f} eV
"""
    print(result)

    log_path = os.path.join(args.aexx_dir, 'aexx_fit.log')
    with open(log_path, 'w') as f:
        f.write(result)
    print(f"Results saved to: {log_path}")

if __name__ == '__main__':
    main()

Usage:

python fit_aexx.py --exp_gap 8.8 --aexx_dir Bulk/AEXX/

Step D: Apply Fitted AEXX

After fitting, update the AEXX value in all HSE INCAR templates. The fitted value replaces the default 0.25 in:

• write_incar_relax_hse() (LEVEL3 only)
• write_incar_static_hse() (LEVEL2/3)

Log the decision:

[YYYY-MM-DD HH:MM:SS] STEP: AEXX Determination
  Description: AEXX fitting complete
  Method: Linear fit to 3 data points
  Experimental gap: 8.8 eV (source: DOI:xxx)
  Fitted AEXX: 0.32
  Used in: all HSE INCAR files

IMPORTANT: The AEXX fitting (Bulk/AEXX/) can run in parallel with the supercell Bulk/relax/ step. It uses the primitive cell, so it has no dependency on the supercell relaxation. The fitted AEXX value is only needed before writing the first HSE INCAR file (statics for LEVEL2, or relax for LEVEL3).

Workflow Integration

For LEVEL2:
  Bulk/relax (PBE) ──────────────────────────→ defect relax (PBE) → statics (HSE, AEXX fitted)
  Bulk/ktest → Bulk/AEXX (parallel, prim) ──→ fit AEXX ──────────────↗

For LEVEL3:
  Bulk/ktest → Bulk/AEXX (prim) → fit AEXX → Bulk/relax (HSE) → defect relax (HSE) → statics (HSE)

For LEVEL3, the AEXX fitting must complete before Bulk/relax starts, since Bulk/relax itself uses HSE.

K-mesh Determination

When determining k-mesh, you MUST ask the user to choose one of the following options:

1. Automatic (recommended) — Auto-determine k-mesh using k-spacing formula based on cell size
2. User-specified — User provides a specific k-mesh (e.g., 4×4×2)
3. Convergence test — Run calculations at multiple k-mesh densities to verify convergence

CRITICAL: For any calculations on small cells (primitive cells, AEXX fitting, etc.), convergence testing is MANDATORY. Only supercell calculations (60+ atoms) can safely use Gamma-only without testing.

N×a ≈ R Method (Auto Mode)

k-mesh density is determined by the parameter R (in Å): for each lattice direction, N_i = max(1, round(R / a_i)), where a_i is the lattice constant along that direction. Larger R = denser mesh.

Recommended R values:

Example (Al₂O₃, a=4.76 Å, c=12.99 Å):

Python helper function:

import numpy as np

def get_kmesh(lattice_vectors, R=30):
    """Determine k-mesh from lattice vectors and target R parameter.
    lattice_vectors: 3x3 array (rows are vectors, in Å)
    R: target N*a product in Å (larger = denser)
    Returns: [N1, N2, N3]
    """
    kmesh = []
    for i in range(3):
        a_len = np.linalg.norm(lattice_vectors[i])
        n = max(1, round(R / a_len))
        kmesh.append(n)
    return kmesh

Supercell K-mesh (Automatic — No Test Needed)

For supercells with 60+ atoms (lattice vectors > ~10 Å), Gamma-only (1×1×1) is almost always sufficient:

Gamma
0
G
1 1 1
0 0 0

Rationale: Reciprocal vectors are very short (e.g., 2π/13 ≈ 0.48 Å⁻¹), so a single k-point adequately samples the BZ. Use vasp_gam for speed.

Primitive Cell K-mesh (Convergence Test Required)

For AEXX fitting, band structure, DOS, or any property calculated on the primitive cell, a k-mesh convergence test is MANDATORY.

Convergence test procedure:

1. Prepare test calculations in Bulk/ktest/ with increasing k-mesh density:

Bulk/ktest/
├── k1/    # e.g., 3×3×1 (coarse, R≈15)
│   ├── POSCAR INCAR KPOINTS POTCAR submit_job
├── k2/    # e.g., 4×4×2 (medium, R≈20)
│   ├── ...
├── k3/    # e.g., 6×6×2 (dense, R≈30)
│   ├── ...
├── k4/    # e.g., 8×8×3 (very dense, R≈40, optional)
│   ├── ...
└── ktest_result.log

2. Choose test k-meshes based on the lattice symmetry and cell dimensions:

   • Start from R≈15 (coarse), increase to R≈20, 30, 40
   • Ensure meshes respect lattice symmetry (e.g., hexagonal: N1=N2≠N3)
   • For HSE: fewer points needed (HSE converges faster with k-mesh), R≈20–30 usually sufficient

3. Run PBE static at each k-mesh (fast, even if final target is HSE). Compare total energy.

4. Convergence criterion: Total energy change < 1 meV/atom between successive meshes.

5. Select the converged k-mesh and use it for all subsequent primitive cell calculations (AEXX fitting, etc.).

ktest INCAR (PBE static on primitive cell):

NELM   = 100
EDIFF  = 1E-6
PREC   = Accurate
ISMEAR = 0
SIGMA  = 0.05
GGA    = PE
IALGO  = 38
LORBIT = 10
ISPIN  = 2
ENCUT  = <same as supercell>
NPAR   = 4
LWAVE  = .FALSE.
LCHARG = .FALSE.

Log example:

[YYYY-MM-DD HH:MM:SS] STEP: K-mesh Convergence Test
  Description: 对 {material} 原胞进行 k-mesh 收敛测试
  Cell: a=4.76 Å, b=4.76 Å, c=12.99 Å (rhombohedral, 10 atoms)
  Test results:
    2×2×1: E = -XXX.XXXX eV  (ΔE = --- meV/atom)
    4×4×2: E = -XXX.XXXX eV  (ΔE = X.X meV/atom)
    6×6×3: E = -XXX.XXXX eV  (ΔE = 0.X meV/atom)  ← CONVERGED
    8×8×4: E = -XXX.XXXX eV  (ΔE = 0.0X meV/atom)
  Selected k-mesh: 6×6×3
  Reasoning: Energy converged to < 1 meV/atom

Workflow Integration (K-mesh Test + AEXX Fitting)

For LEVEL2/3 with AEXX fitting (Option 3), the k-mesh convergence test should run before the AEXX fitting calculations:

Bulk/ktest/ (PBE static, fast) → determine converged k-mesh
    ↓
Bulk/AEXX/aexx_*/  (HSE static, use converged k-mesh) → fit AEXX

Both can run in parallel with the supercell workflow (Bulk/relax, defect relax) since they use the primitive cell.

Workflow Overview

CRITICAL: 各步骤存在严格的前后依赖关系，不能提前准备后续步骤的输入文件。

Step 0 (LEVEL2/3): K-mesh Test + AEXX Determination
    询问用户 k-mesh: 自动 / 手动指定 / 收敛测试
    询问用户 AEXX: 标准HSE06(0.25) / 自定义AEXX / 拟合实验带隙
    若选拟合: k-mesh收敛测试(Bulk/ktest/)→搜索文献→AEXX拟合(Bulk/AEXX/)→确定AEXX
    LEVEL2: 可与 Step 1-4 并行（AEXX 在 statics 前确定即可）
    LEVEL3: 必须在 Step 2 前完成（relax 已使用 HSE）
    ↓
Step 1: Supercell Construction
    primitive POSCAR → nearly-cubic supercell (min_atom ~ max_atom)
    → 只准备 Bulk/relax/ 的输入文件 (POSCAR, INCAR, KPOINTS, POTCAR, submit_job)
    ↓
Step 2: Bulk Relaxation (ISIF=3)
    LEVEL1/2: PBE relax | LEVEL3: HSE relax (with fitted AEXX)
    提交 Bulk/relax/ 任务，等待收敛
    ↓ 【必须等 Bulk/relax 完成】
Step 3: Defect Structure Generation + Neutral Relaxation Setup
    用 Bulk/relax/CONTCAR 构建各缺陷 POSCAR
    → 准备所有 Defect/*/q0/relax/ 的输入文件
    ↓
Step 4: Neutral Defect Relaxation (ISIF=2, no NELECT)
    LEVEL1/2: PBE relax | LEVEL3: HSE relax
    提交所有 q0/relax/ 任务（可并行），等待收敛
    ↓ 【必须等所有 q0/relax 完成】
Step 5: Preliminary Statics (Bulk + q0 only) + Charge State Determination
    LEVEL1: PBE statics | LEVEL2/3: HSE statics
    先仅对 Bulk 和所有 q0 进行 statics 计算
    → Bulk/statics/ + 所有 Defect/*/q0/statics/
    ↓ 【必须等 Bulk/statics 和 q0/statics 完成】
    分析 q0/statics/EIGENVAL 与 Bulk/statics/EIGENVAL 对比确定各缺陷的带电态
    ⚠ 必须使用 statics/EIGENVAL（而非 relax/EIGENVAL）：
      LEVEL1: statics 与 relax 均为 PBE，差异不大，但 statics 更准确
      LEVEL2: statics 为 HSE，带隙和缺陷能级位置与 PBE relax 显著不同
      LEVEL3: statics 为 HSE，与 HSE relax 一致，但 statics 更严格
    ↓
Step 6: Charged Defect Relaxation Setup + Submission
    LEVEL1/2: PBE relax | LEVEL3: HSE relax
    用 q0/relax/CONTCAR 作为起始结构
    → 准备所有 q±n/relax/ 的输入文件，提交任务
    ↓ 【必须等所有 charged relax 完成】
Step 7: Charged Defect Static Calculations
    LEVEL1: PBE statics | LEVEL2/3: HSE statics
    用各 charged relax/CONTCAR 准备 statics/ 输入文件
    → 所有 Defect/*/q±n/statics/
    （Bulk/statics 和 q0/statics 在 Step 5 已完成，无需重复）
    ↓ 【必须等所有 statics 完成】
Step 8: Finite-Size Corrections
    FNV (Freysoldt) or Kumagai correction for charged defects
    ↓
Step 9: Formation Energy & Diagram
    compute E_f, plot formation energy vs Fermi level

输入文件准备时机与结构来源

关键：带电缺陷的起始结构是中性弛豫后的 CONTCAR，不是从 bulk 重新构建的缺陷结构。 这样做的原因是中性缺陷弛豫后原子位置已经松弛到缺陷附近的平衡位置，以此为起点加/减电子再弛豫，收敛更快、结果更可靠。

Step 1: Supercell Construction

Goal: Build a nearly-cubic supercell with 60–70 atoms (configurable via min_atom/max_atom).

Source: 原胞结构的获取途径（必须记录到日志）：

1. 用户直接提供 — 已有的 POSCAR 或 CIF 文件
2. 数据库下载 — Materials Project、ICSD、COD 等，需记录材料ID（如 mp-406）
3. 文献构建 — 根据实验晶格参数手动构建，需记录文献来源
4. 其他来源 — 如从已有计算结果中提取

Method:

1. Read the primitive cell lattice vectors
2. Find a transformation matrix P such that P × primitive_cell ≈ cubic, with N_atoms in [min_atom, max_atom]
3. The supercell should minimize the ratio (max_lattice_length / min_lattice_length) to be as cubic as possible
4. For CdTe zincblende (FCC primitive, 2 atoms): a 4×4×4 FCC = 2×2×2 conventional cubic = 64 atoms

Log: Record:

• 原胞来源： 从哪获取（数据库名+材料ID / 文献DOI / 用户提供的文件路径）
• 原胞信息： 空间群、晶格参数、原子数
• 超胞变换矩阵、目标原子数、实际原子数
• 超胞晶格向量、cubicity ratio

Example supercell (CdTe 64 atoms):

Cubic_cell
1.0
13.258  0.000  0.000
 0.000 13.258  0.000
 0.000  0.000 13.258
Cd Te
32 32

Step 2: Bulk Relaxation (ISIF=3)

此步骤是整个流程的起点。初始阶段只准备 Bulk/relax/ 目录的输入文件。

输入文件准备： 在 Bulk/relax/ 中放置 POSCAR, INCAR, KPOINTS, POTCAR, submit_job

INCAR for bulk relaxation (PBE — LEVEL1/2):

NSW = 500
NELM = 100
ISIF = 3              # Full cell + ion relaxation
IBRION = 2
POTIM = 0.5
EDIFF = 1E-5
EDIFFG = -0.01
ISMEAR = 0
SIGMA = 0.05
GGA = PE              # PBE functional
IALGO = 38
LORBIT = 10
PREC = Accurate
ISPIN = 2
ENCUT = <1.3×ENMAX>   # From POTCAR
NPAR = 8

INCAR for bulk relaxation (HSE06 — LEVEL3):

NSW = 500
NELM = 100
ISIF = 3              # Full cell + ion relaxation
IBRION = 2
POTIM = 0.3           # Smaller for HSE stability
EDIFF = 1E-5
EDIFFG = -0.02        # Slightly relaxed for HSE
ISMEAR = 0
SIGMA = 0.05
GGA = PE
LHFCALC  = .TRUE.     # Enable Hartree-Fock
HFSCREEN = 0.2         # HSE06 screening
AEXX     = 0.25        # 25% exact exchange
ALGO     = Damped      # Required for HF relaxation (NOT IALGO=38)
TIME     = 0.4
PRECFOCK = Fast
LREAL  = .FALSE.       # Required for HF accuracy
LORBIT = 10
PREC = Accurate
ISPIN = 2
ENCUT = <1.3×ENMAX>
NPAR = 8

KPOINTS: Gamma-only for large supercells (64+ atoms):

Gamma
0
G
1 1 1
0 0 0

After relaxation:

• Use CONTCAR as the relaxed bulk structure
• Record the final energy from OSZICAR
• This CONTCAR is the basis for ALL defect structures

Log: Record ENCUT source, k-mesh, final energy, lattice parameters change, convergence status.

Step 3: Defect Structure Generation（Bulk弛豫完成后）

前置条件： Bulk/relax/ 已收敛，CONTCAR 可用。

操作： 用 Bulk/relax/CONTCAR 构建各缺陷 POSCAR，并准备所有 q0/relax/ 目录的完整输入文件（POSCAR, INCAR, KPOINTS, POTCAR, submit_job）。

Starting from the relaxed bulk CONTCAR, create defect POSCARs:

Vacancies

Remove one atom from the supercell:

• V_Cd: Remove 1 Cd atom → (N_Cd - 1) Cd + N_Te Te
• V_Te: Remove 1 Te atom → N_Cd Cd + (N_Te - 1) Te

Antisites

Replace one atom with another species:

• Te_Cd: Replace 1 Cd with Te → (N_Cd - 1) Cd + (N_Te + 1) Te
• Cd_Te: Replace 1 Te with Cd → (N_Cd + 1) Cd + (N_Te - 1) Te

Impurities (Substitutional)

Replace one host atom with an impurity:

• Cl_Te: Replace 1 Te with Cl → N_Cd Cd + (N_Te - 1) Te + 1 Cl
• Na_Cd: Replace 1 Cd with Na → (N_Cd - 1) Cd + N_Te Te + 1 Na

Defect complexes

Create two point defects in the same supercell at various nearest-neighbor distances (nn1, nn2, nn3).

Log for EACH defect:

• Which atom was removed/replaced (index, species, fractional coordinates)
• Resulting composition and atom count
• Defect position in the supercell

Step 4: Neutral Defect Relaxation (ISIF=2)

输入文件已在 Step 3 中准备好（在 q0/relax/ 目录下）。

INCAR for neutral defect relaxation (PBE — LEVEL1/2):

NSW = 500
NELM = 100
ISIF = 2              # Fix cell, relax ions only
IBRION = 2
POTIM = 0.5
EDIFF = 1E-5
EDIFFG = -0.01
ISMEAR = 0
SIGMA = 0.05
GGA = PE
IALGO = 38
LORBIT = 10
PREC = Accurate
ISPIN = 2
ENCUT = <same as bulk>
NPAR = 8

INCAR for neutral defect relaxation (HSE06 — LEVEL3):

NSW = 500
NELM = 100
ISIF = 2              # Fix cell, relax ions only
IBRION = 2
POTIM = 0.3
EDIFF = 1E-5
EDIFFG = -0.02
ISMEAR = 0
SIGMA = 0.05
GGA = PE
LHFCALC  = .TRUE.
HFSCREEN = 0.2
AEXX     = 0.25
ALGO     = Damped
TIME     = 0.4
PRECFOCK = Fast
LREAL  = .FALSE.
LORBIT = 10
PREC = Accurate
ISPIN = 2
ENCUT = <same as bulk>
NPAR = 8

Key: No NELECT tag → VASP auto-determines from POTCAR (neutral defect).

After relaxation: CONTCAR becomes the starting structure for charged states.

Log: Final energy, convergence status, max force, number of ionic steps.

Step 5: Preliminary Statics + Charge State Determination

前置条件： Bulk/relax/ 和所有 q0/relax/ 已收敛。

操作顺序：

1. 先对 Bulk 和所有 q0 缺陷运行 statics 计算（LEVEL1: PBE, LEVEL2/3: HSE）
2. 等待 statics 收敛后，用 statics/EIGENVAL 进行电荷态分析

⚠ 关键：必须使用 statics/EIGENVAL 而非 relax/EIGENVAL：

• relax/EIGENVAL 来自弛豫过程中的电子结构，精度较低（NSW≠0，离子仍在移动）
• statics/EIGENVAL 来自固定结构的精确单点计算，VBM/CBM/缺陷能级位置更准确
• LEVEL2 尤其关键：relax 用 PBE（带隙偏小），statics 用 HSE（带隙接近实验值），两者的缺陷能级位置可能有本质差异，直接影响电荷态判断

Core Method: Compare Defect vs Bulk EIGENVAL

在 Bulk/statics/ 和 q0/statics/ 完成后，比较它们的 EIGENVAL 文件：

1. Identify the bulk band structure特征 — VBM, CBM, 以及价带/导带的能量分布模式
2. 在缺陷EIGENVAL中找到缺陷能级 — 不属于bulk价带或导带的"多出来"的态
3. 分析缺陷能级的占据情况
4. 据此确定可能的带电态

Step 5a: 建立bulk参考 — 提取带边和能带分布特征

EIGENVAL format (ISPIN=2): band_index E_up E_down occ_up occ_down

首先从host EIGENVAL读取带边附近的完整能带结构：

Host (CdTe 64-atom, NELECT=576, 288 e⁻/spin):

  --- 价带 (全部 occ=1.0) ---
  Band 276-277: E ≈ 1.081 eV   ← 价带态 (2重简并)
  Band 278-285: E ≈ 1.108 eV   ← 价带态 (8重简并)
  Band 286-288: E ≈ 1.790 eV   ← VBM (3重简并，价带顶)
  ─────────── 带隙 0.57 eV ───────────
  Band 289:     E ≈ 2.358 eV   ← CBM (导带底)
  Band 290-293: E ≈ 3.535 eV   ← 导带态 (4重简并)
  --- 导带 (全部 occ=0.0) ---

关键特征要记住：

• 价带顶(VBM): 1.790 eV (3重简并)
• 导带底(CBM): 2.358 eV
• 带隙: 0.57 eV
• 价带中的简并模式和能量间距

Step 5b: 在缺陷EIGENVAL中识别缺陷能级

识别方法 — 三个判据（综合使用）：

判据1: 部分占据 (最直接) 占据数在0和1之间的态 (如 occ=0.667) 一定是缺陷能级。因为bulk中所有态要么完全占据(1.0)要么完全空(0.0)，部分占据只会出现在费米面附近的缺陷态上。

判据2: 在bulk带隙内出现的态 能量在bulk VBM (1.790) 和 CBM (2.358) 之间的态是缺陷能级。但要注意：超胞中的本体态会因缺陷影响略有偏移（~0.1 eV），所以需要一定容差。

判据3: 偏离bulk能带分布模式的"额外"态 对比bulk和缺陷的能带分布，找到能量位置明显不同于任何bulk能带的态。具体做法：

• bulk中价带顶是3重简并的1.790 eV，下面是8重简并的1.108 eV
• 如果在缺陷EIGENVAL中，在1.108和1.790之间出现了额外的态，或在2.358上方出现了被占据的态，这些就是缺陷态

Step 5b 实例详解

V_Cd 中性 (NELECT=564, 282 e⁻/spin):

  Band 275-277: E ≈ 1.128 eV, occ=1.0  ← 对应bulk的 ~1.108价带态 (轻微偏移)
  Band 278-280: E ≈ 1.360 eV, occ=1.0  ← 对应bulk的 ~1.790 VBM? 不对!
                                           bulk VBM是3重简并的1.790
                                           这里变成了1.360, 差了0.43 eV
                                           → 这是价带态,只是被缺陷推移了
  Band 281-283: E ≈ 1.751 eV, occ=0.667 ←←← 缺陷能级!!!
                                           ✓ 判据1: 部分占据(0.667≠0或1)
                                           ✓ 判据2: 在bulk gap附近(1.75, 接近VBM 1.79)
                                           ✓ 判据3: bulk中此位置没有部分占据的3重简并态
  Band 284:     E ≈ 2.382 eV, occ=0.0   ← 对应bulk CBM(2.358), 导带态

V_Te 中性 (NELECT=570, 285 e⁻/spin):

  Band 282-284: E ≈ 1.593 eV, occ=1.0  ← 对应bulk ~1.790 VBM (下移了0.2 eV,正常范围)
  Band 285:     E ≈ 2.187 eV, occ=1.0  ←←← 缺陷能级!!!
                                           ✓ 判据2: 在bulk带隙内 (1.790 < 2.187 < 2.358)
                                           ✓ 判据3: bulk中VBM(3重简并)和CBM之间无态
                                              而这里多出了一个完全占据的态
                                           注意: occ=1.0, 但位置暴露了它
  Band 286:     E ≈ 3.320 eV, occ=0.0  ← 导带态

Te_Cd 中性 (NELECT=570, 285 e⁻/spin):

  Band 282-284: E ≈ 1.849 eV, occ=1.0  ← 对应bulk VBM (~1.790, 轻微上移)
  Band 285:     E ≈ 2.566 eV, occ≈0.82 ←←← 缺陷能级!!!
                                           ✓ 判据1: 部分占据
                                           ✓ 判据3: 在bulk CBM(2.358)上方, 却有电子!
  Band 286-288: E ≈ 2.64-2.72 eV, occ≈0.1 ←←← 缺陷能级!!!
                                           ✓ 判据1: 部分占据
                                           ✓ 判据3: 比CBM还高, 却有残余电子
  Band 289:     E ≈ 3.698 eV, occ=0.0  ← 导带态

Cd_Te 中性 (NELECT=582, 291 e⁻/spin):

  Band 288-290: E ≈ 1.643 eV, occ=1.0  ← 对应bulk VBM
  Band 291:     E ≈ 2.358 eV, occ=1.0  ←←← 缺陷能级!!!
                                           ✓ 判据2: 恰好在bulk CBM能量(2.358 eV)
                                           ✓ 判据3: bulk中此处是空的CBM,
                                              现在却完全占据→额外的电子在此
  Band 292:     E ≈ 3.327 eV, occ=0.0  ← 导带态

Step 5c: 判断方法总结 — 实操流程

1. 从host EIGENVAL提取:
   - VBM能量和简并度
   - CBM能量
   - 价带和导带的能量分布模式(各简并组的能量位置)

2. 在defect EIGENVAL中扫描:
   (a) 先找部分占据态 → 一定是缺陷能级 [判据1, 最可靠]
   (b) 再检查带隙区间(VBM~CBM)内有无额外态 → 缺陷能级 [判据2]
   (c) 对比能带分布模式,找不属于价带或导带的态 → 缺陷能级 [判据3]

3. 统计缺陷能级中的电子数:
   - 总电子数 = Σ (occ_up + occ_down) for all defect bands

4. 确定电荷态:
   - 缺陷能级中有N_occ个电子, N_empty个空态
   - 可移走的电子 → 正电荷态: q = +1, +2, ..., +N_occ
   - 可加入的电子 → 负电荷态: q = -1, -2, ..., -N_empty

Step 5d: 从缺陷能级中的电子数确定电荷态和NELECT

General Rules for Charge State Determination

1. Identify defect levels: States inside or near the bulk gap that differ from the bulk band structure
2. Count electrons in defect levels:
   • Each band holds 1 electron per spin (2 total for ISPIN=2)
   • Partial occupation (e.g., 0.667) means fractional filling
3. Determine removable electrons: Electrons in defect states above VBM can be removed → positive q
4. Determine addable electrons: Empty defect states below CBM can accept electrons → negative q
5. Set NELECT:

   NELECT = NELECT_neutral - q
   

   (q>0 means fewer electrons; q<0 means more electrons)

NELECT Reference Table (CdTe)

Log for EACH defect: Record:

• Bulk VBM and CBM eigenvalues
• Defect level positions and occupations (band index, energy, occ_up, occ_down)
• Number of electrons in defect levels
• Reasoning: which electrons can be added/removed
• Resulting charge state list and NELECT values

Step 6: Charged Defect Relaxation（q0弛豫完成后）

前置条件： 所有 q0/relax/ 已收敛，且 Step 5 的电荷态分析已完成。

操作：

1. 将 q0/relax/CONTCAR 复制为各 q±n/relax/POSCAR
2. 准备 INCAR（加 NELECT）、KPOINTS、POTCAR、submit_job
3. 提交所有 charged relax 任务（可并行）

Starting structure: CONTCAR from the neutral (q=0) relaxation.

INCAR: Same as neutral relaxation, but add:

NELECT = <calculated value>

Directory structure:

ProjectRoot/
├── Bulk/
│   ├── ktest/                    # K-mesh convergence test (primitive cell)
│   │   ├── primitive_POSCAR      # Primitive cell
│   │   ├── k1/                   # Coarse k-mesh (e.g., 2×2×1)
│   │   ├── k2/                   # Medium k-mesh (e.g., 4×4×2)
│   │   ├── k3/                   # Dense k-mesh (e.g., 6×6×3)
│   │   └── ktest_result.log      # Convergence results
│   ├── AEXX/                     # AEXX fitting (LEVEL2/3 only)
│   │   ├── primitive_POSCAR      # Primitive cell for fitting
│   │   ├── aexx_0.20/            # HSE static with AEXX=0.20
│   │   ├── aexx_0.25/            # HSE static with AEXX=0.25
│   │   ├── aexx_0.30/            # HSE static with AEXX=0.30
│   │   ├── fit_aexx.py           # Fitting script
│   │   └── aexx_fit.log          # Fitting results
│   ├── relax/                    # Bulk relaxation (ISIF=3)
│   └── statics/                  # Bulk static calculation
├── Defect/
│   ├── Intrinsic/
│   │   ├── V_Cd/
│   │   │   ├── q0/
│   │   │   │   ├── relax/        # Neutral relaxation
│   │   │   │   └── statics/      # Static after relax
│   │   │   ├── q-1/
│   │   │   │   ├── relax/        # Charged relaxation (NELECT+1)
│   │   │   │   └── statics/
│   │   │   └── q-2/
│   │   │       ├── relax/        # Charged relaxation (NELECT+2)
│   │   │       └── statics/
│   │   ├── V_Te/
│   │   │   ├── q0/
│   │   │   │   ├── relax/
│   │   │   │   └── statics/
│   │   │   ├── q+1/
│   │   │   │   ├── relax/
│   │   │   │   └── statics/
│   │   │   └── q+2/
│   │   │       ├── relax/
│   │   │       └── statics/
│   │   ├── Te_Cd/
│   │   │   └── ...             # Same pattern: q*/relax/, q*/statics/
│   │   └── Cd_Te/
│   │       └── ...
│   └── Impurity/
│       ├── Cl_Te/
│       │   └── ...             # Same pattern
│       └── Na_Cd/
│           └── ...
├── ThermoStability/              # Chemical potential & stability analysis
└── FormationEnergy/              # Formation energy results & plots

Key: Each charged calculation starts from the relaxed neutral CONTCAR, NOT from the unrelaxed defect POSCAR.

Log: Record starting POSCAR source, NELECT value, convergence for each charge state.

Step 7: Static Calculations（所有弛豫完成后）

前置条件： Bulk/relax/ 和所有 Defect/*/q*/relax/ 已收敛。

操作：

1. 将各 relax/CONTCAR 复制为对应 statics/POSCAR
2. 将各 relax/CHGCAR 复制到对应 statics/（ICHARG=1 需要读取）
3. 准备 INCAR（静态参数）、KPOINTS、POTCAR、submit_job
4. 提交所有 statics 任务（可并行）

需要准备的 statics 目录：

• Bulk/statics/ — 使用 Bulk/relax/CONTCAR
• 所有 Defect/*/q*/statics/ — 使用对应 relax/CONTCAR

INCAR for static (PBE — LEVEL1):

NELM = 100
EDIFF = 1E-5
ISMEAR = 0
SIGMA = 0.05
GGA = PE
IALGO = 38
LREAL = Auto
LORBIT = 10
PREC = Accurate
ISPIN = 2
ENCUT = <same as relaxation>
NPAR = 8
ICORELEVEL = 1        # Core-level eigenvalues for Kumagai correction
LVHAR = .TRUE.        # Write LOCPOT for Freysoldt correction
ICHARG = 1            # Read CHGCAR from relaxation

INCAR for static (HSE06 — LEVEL2/3):

NELM = 200            # HSE needs more SCF steps
EDIFF = 1E-5
ISMEAR = 0
SIGMA = 0.05
GGA = PE
LHFCALC  = .TRUE.
HFSCREEN = 0.2
AEXX     = 0.25
ALGO     = All        # Use All for HSE static (NOT Damped)
TIME     = 0.4
PRECFOCK = Fast
LREAL  = .FALSE.      # Required for HF accuracy
LORBIT = 10
PREC = Accurate
ISPIN = 2
ENCUT = <same as relaxation>
NPAR = 8
ICORELEVEL = 1
LVHAR = .TRUE.
ICHARG = 1

Note: No NSW, no IBRION → single-point calculation. Copy CHGCAR from the relaxation to the static directory (ICHARG=1 reads it).

Required outputs: LOCPOT (for FNV correction), OUTCAR (for energies and site potentials).

Log: Record total energies from each static calculation.

Step 8: Finite-Size Corrections（所有statics完成后）

前置条件： Bulk/statics/ 和所有 Defect/*/q*/statics/ 已收敛。

fnv_correction.py — FNV 修正计算脚本

#!/usr/bin/env python3
"""
Freysoldt (FNV) correction for charged defects.
Usage: python fnv_correction.py <bulk_locpot> <defect_locpot> <charge> <epsilon> [--madelung <alpha>]
Output: correction value in eV, writes detailed log
"""
import numpy as np
import sys
import os
import argparse

def read_locpot(filepath):
    """读取VASP LOCPOT文件，返回晶格向量和3D势场数据"""
    with open(filepath) as f:
        f.readline()  # comment
        scale = float(f.readline())
        lattice = np.array([[float(x) for x in f.readline().split()] for _ in range(3)]) * scale
        species = f.readline().split()
        counts = [int(x) for x in f.readline().split()]
        total_atoms = sum(counts)
        # Skip coordinate type line and atom coordinates
        f.readline()  # Direct/Cartesian
        for _ in range(total_atoms):
            f.readline()
        f.readline()  # blank line
        # Read grid dimensions
        nx, ny, nz = [int(x) for x in f.readline().split()]
        # Read potential data
        data = []
        while len(data) < nx * ny * nz:
            line = f.readline()
            if line.strip():
                data.extend([float(x) for x in line.split()])
        pot = np.array(data).reshape((nz, ny, nx)).transpose((2, 1, 0))  # [nx,ny,nz]
    return lattice, pot, (nx, ny, nz)

def planar_average(pot, axis):
    """沿指定轴计算平面平均势"""
    axes = [0, 1, 2]
    axes.remove(axis)
    return np.mean(pot, axis=tuple(axes))

def compute_madelung_energy(alpha_M, q, epsilon, L):
    """计算 Madelung 能量 (eV)
    E_Mad = alpha_M * q^2 / (2 * epsilon * L)
    L in Angstrom, result in eV (with conversion factor)
    """
    # Hartree atomic units: 1 Ha = 27.2114 eV, 1 Bohr = 0.529177 Å
    L_bohr = L / 0.529177
    E_hartree = alpha_M * q**2 / (2.0 * epsilon * L_bohr)
    return E_hartree * 27.2114

def compute_potential_alignment(bulk_pot, defect_pot, lattice, axis, q, epsilon):
    """计算沿某轴的势能对齐量 ΔV"""
    V_bulk = planar_average(bulk_pot, axis)
    V_defect = planar_average(defect_pot, axis)
    diff = V_defect - V_bulk

    n = len(diff)
    # 取远离缺陷的区域（两端各1/4）作为对齐参考
    far_region = np.concatenate([diff[:n//4], diff[3*n//4:]])
    delta_V = np.mean(far_region)
    return delta_V

def main():
    parser = argparse.ArgumentParser(description='FNV correction for charged defects')
    parser.add_argument('bulk_locpot', help='Path to bulk LOCPOT')
    parser.add_argument('defect_locpot', help='Path to defect LOCPOT')
    parser.add_argument('charge', type=int, help='Charge state q')
    parser.add_argument('epsilon', type=float, help='Dielectric constant')
    parser.add_argument('--madelung', type=float, default=2.8373,
                        help='Madelung constant (default: 2.8373 for SC)')
    parser.add_argument('--output', default=None, help='Output log file')
    args = parser.parse_args()

    q = args.charge
    if q == 0:
        print("q=0: No correction needed")
        print("CORRECTION: 0.000")
        return

    print(f"Reading bulk LOCPOT: {args.bulk_locpot}")
    lattice_b, pot_b, grid_b = read_locpot(args.bulk_locpot)
    print(f"Reading defect LOCPOT: {args.defect_locpot}")
    lattice_d, pot_d, grid_d = read_locpot(args.defect_locpot)

    # 超胞特征长度（取晶格向量长度的平均值）
    L_values = [np.linalg.norm(lattice_b[i]) for i in range(3)]
    L_avg = np.mean(L_values)

    # Madelung energy
    E_mad = compute_madelung_energy(args.madelung, q, args.epsilon, L_avg)

    # Potential alignment（三个轴取平均）
    delta_Vs = []
    for axis in range(3):
        dV = compute_potential_alignment(pot_b, pot_d, lattice_b, axis, q, args.epsilon)
        delta_Vs.append(dV)
    delta_V = np.mean(delta_Vs)

    # Total correction
    E_align = q * delta_V
    E_corr = E_mad + E_align

    # Output
    result = f"""FNV Correction Result:
  Charge state: q = {q}
  Dielectric constant: epsilon = {args.epsilon}
  Madelung constant: alpha = {args.madelung}
  Supercell lengths: {L_values[0]:.3f}, {L_values[1]:.3f}, {L_values[2]:.3f} Å
  Average length: {L_avg:.3f} Å
  ---
  Madelung energy: {E_mad:.4f} eV
  Potential alignment (per axis): {delta_Vs[0]:.4f}, {delta_Vs[1]:.4f}, {delta_Vs[2]:.4f} eV
  Average ΔV: {delta_V:.4f} eV
  Alignment correction (q×ΔV): {E_align:.4f} eV
  ---
  Total FNV correction: {E_corr:.4f} eV
CORRECTION: {E_corr:.4f}"""

    print(result)
    if args.output:
        with open(args.output, 'w') as f:
            f.write(result + '\n')

if __name__ == '__main__':
    main()

使用方法

# 对每个带电缺陷运行
python fnv_correction.py \
    Bulk/statics/LOCPOT \
    Defect/Intrinsic/V_Cd/q-2/statics/LOCPOT \
    -2 \
    10.3 \
    --madelung 2.8373 \
    --output Defect/Intrinsic/V_Cd/q-2/fnv_correction.log

Step 9: Formation Energy & Diagram

formation_energy.py — 形成能计算与绘图脚本

#!/usr/bin/env python3
"""
Formation energy calculation and diagram plotting.
Reads VASP results, applies corrections, plots E_f vs E_Fermi.

Usage:
  python formation_energy.py --config formation_config.yaml
  python formation_energy.py --workdir . --epsilon 10.3 --band_gap 1.5

Config file (formation_config.yaml):
  bulk_static: Bulk/statics
  epsilon: 10.3
  madelung: 2.8373
  band_gap: 1.5
  E_pure:
    Cd: -1.7736
    Te: -4.6974
  chemical_potential_limits:
    Cd-rich:
      Cd: 0.0
      Te: -1.1854
    Te-rich:
      Cd: -1.1854
      Te: 0.0
  defects:
    - name: V_Cd
      path: Defect/Intrinsic/V_Cd
      removed: {Cd: 1}
      added: {}
    - name: V_Te
      path: Defect/Intrinsic/V_Te
      removed: {Te: 1}
      added: {}
    - name: Te_Cd
      path: Defect/Intrinsic/Te_Cd
      removed: {Cd: 1}
      added: {Te: 1}
    - name: Cl_Te
      path: Defect/Impurity/Cl_Te
      removed: {Te: 1}
      added: {Cl: 1}
"""
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
import yaml
import subprocess
import re
from datetime import datetime

def get_energy_from_outcar(outcar_path):
    """从OUTCAR提取 'energy without entropy' (sigma->0)"""
    E = None
    with open(outcar_path) as f:
        for line in f:
            if "energy  without entropy" in line:
                E = float(line.split()[-1])
    return E

def get_vbm_from_eigenval(eigenval_path):
    """从EIGENVAL提取VBM（最高占据态能量）"""
    with open(eigenval_path) as f:
        lines = f.readlines()
    # Header: line 6 has NELECT, NKPTS, NBANDS
    header = lines[5].split()
    nelect = int(header[0])
    nkpts = int(header[1])
    nbands = int(header[2])

    vbm = -999.0
    idx = 7  # skip header lines
    for k in range(nkpts):
        idx += 1  # blank line
        idx += 1  # k-point header
        for b in range(nbands):
            parts = lines[idx].split()
            band_idx = int(parts[0])
            e_up = float(parts[1])
            occ_up = float(parts[2])
            if len(parts) > 3:  # ISPIN=2
                e_down = float(parts[3])
                occ_down = float(parts[4])
                if occ_up > 0.5:
                    vbm = max(vbm, e_up)
                if occ_down > 0.5:
                    vbm = max(vbm, e_down)
            else:
                if occ_up > 0.5:
                    vbm = max(vbm, e_up)
            idx += 1
    return vbm

def get_nelect_from_outcar(outcar_path):
    """从OUTCAR提取NELECT"""
    with open(outcar_path) as f:
        for line in f:
            if "NELECT" in line:
                return int(float(line.split()[2]))
    return None

def compute_fnv_correction(bulk_locpot, defect_locpot, charge, epsilon, madelung=2.8373):
    """调用fnv_correction.py并返回修正值"""
    if charge == 0:
        return 0.0
    script = os.path.join(os.path.dirname(__file__), 'fnv_correction.py')
    try:
        result = subprocess.run(
            ['python3', script, bulk_locpot, defect_locpot, str(charge), str(epsilon),
             '--madelung', str(madelung)],
            capture_output=True, text=True, timeout=300
        )
        for line in result.stdout.split('\n'):
            if line.startswith('CORRECTION:'):
                return float(line.split(':')[1].strip())
    except Exception as e:
        print(f"  WARNING: FNV correction failed: {e}")
    return 0.0

def find_charge_dirs(defect_path):
    """查找缺陷目录下的所有电荷态目录"""
    charges = []
    for d in sorted(os.listdir(defect_path)):
        if d.startswith('q') and os.path.isdir(os.path.join(defect_path, d)):
            q_str = d[1:]  # remove 'q'
            try:
                q = int(q_str.replace('+', ''))
                charges.append((q, os.path.join(defect_path, d)))
            except ValueError:
                pass
    return charges

def main():
    import argparse
    parser = argparse.ArgumentParser(description='Formation energy calculation and plotting')
    parser.add_argument('--config', required=True, help='YAML config file')
    parser.add_argument('--output', default='FormationEnergy', help='Output directory')
    args = parser.parse_args()

    with open(args.config) as f:
        cfg = yaml.safe_load(f)

    os.makedirs(args.output, exist_ok=True)
    log_path = os.path.join(args.output, 'formation_energy.log')
    log_lines = []
    def log(msg):
        log_lines.append(f"[{datetime.now().strftime('%Y-%m-%d %H:%M:%S')}] {msg}")
        print(msg)

    # ---- 读取 bulk 能量和 VBM ----
    bulk_dir = cfg['bulk_static']
    E_bulk = get_energy_from_outcar(os.path.join(bulk_dir, 'OUTCAR'))
    E_VBM = get_vbm_from_eigenval(os.path.join(bulk_dir, 'EIGENVAL'))
    nelect_bulk = get_nelect_from_outcar(os.path.join(bulk_dir, 'OUTCAR'))
    band_gap = cfg['band_gap']
    epsilon = cfg['epsilon']
    madelung = cfg.get('madelung', 2.8373)

    log(f"Bulk energy: {E_bulk:.6f} eV")
    log(f"VBM: {E_VBM:.4f} eV")
    log(f"Band gap: {band_gap} eV")
    log(f"Dielectric constant: {epsilon}")

    E_pure = cfg['E_pure']         # {element: E_per_atom}
    limits = cfg['chemical_potential_limits']  # {limit_name: {element: delta_mu}}

    # ---- 计算各缺陷的形成能 ----
    all_results = {}  # {defect_name: {q: {limit: E_f_at_VBM}}}

    for defect_cfg in cfg['defects']:
        name = defect_cfg['name']
        dpath = defect_cfg['path']
        removed = defect_cfg.get('removed', {})
        added = defect_cfg.get('added', {})

        log(f"\n--- Defect: {name} ---")
        all_results[name] = {}

        charge_dirs = find_charge_dirs(dpath)
        for q, q_dir in charge_dirs:
            statics_dir = os.path.join(q_dir, 'statics')
            outcar = os.path.join(statics_dir, 'OUTCAR')
            if not os.path.exists(outcar):
                log(f"  q={q}: OUTCAR not found, skipping")
                continue

            E_defect = get_energy_from_outcar(outcar)
            log(f"  q={q}: E_defect = {E_defect:.6f} eV")

            # FNV correction
            E_corr = 0.0
            if q != 0:
                bulk_locpot = os.path.join(bulk_dir, 'LOCPOT')
                defect_locpot = os.path.join(statics_dir, 'LOCPOT')
                if os.path.exists(bulk_locpot) and os.path.exists(defect_locpot):
                    E_corr = compute_fnv_correction(
                        bulk_locpot, defect_locpot, q, epsilon, madelung)
                    log(f"  q={q}: FNV correction = {E_corr:.4f} eV")
                else:
                    log(f"  q={q}: WARNING - LOCPOT not found, correction=0")

            # 对各化学势条件计算形成能
            all_results[name][q] = {}
            for limit_name, delta_mus in limits.items():
                # 化学势项: Σ n_i × μ_i
                chem_term = 0.0
                for elem, n in removed.items():
                    mu_i = E_pure[elem] + delta_mus[elem]
                    chem_term += n * mu_i
                for elem, n in added.items():
                    mu_i = E_pure[elem] + delta_mus.get(elem, 0.0)
                    chem_term -= n * mu_i

                # E_f at VBM (E_Fermi=0)
                E_f_vbm = E_defect - E_bulk + chem_term + q * E_VBM + E_corr
                all_results[name][q][limit_name] = E_f_vbm
                log(f"  q={q}, {limit_name}: E_f(VBM) = {E_f_vbm:.4f} eV")

    # ---- 计算转变能级 ----
    log(f"\n=== Transition Levels ===")
    for name, q_data in all_results.items():
        charges = sorted(q_data.keys())
        for i in range(len(charges)):
            for j in range(i+1, len(charges)):
                q1, q2 = charges[i], charges[j]
                for limit_name in limits:
                    if limit_name in q_data[q1] and limit_name in q_data[q2]:
                        E1 = q_data[q1][limit_name]
                        E2 = q_data[q2][limit_name]
                        if q2 != q1:
                            E_F_trans = (E1 - E2) / (q2 - q1)
                            if 0 <= E_F_trans <= band_gap:
                                log(f"  {name}: ε({q1}/{q2}) = {E_F_trans:.4f} eV  [{limit_name}]")

    # ---- 绘图 ----
    for limit_name in limits:
        fig, ax = plt.subplots(figsize=(8, 6))
        E_F = np.linspace(0, band_gap, 500)
        colors = plt.cm.tab10(np.linspace(0, 1, len(all_results)))

        for (name, q_data), color in zip(all_results.items(), colors):
            lines = []
            for q, limit_data in sorted(q_data.items()):
                if limit_name in limit_data:
                    E_f_vbm = limit_data[limit_name]
                    E_f_line = E_f_vbm + q * E_F
                    lines.append(E_f_line)
                    # Faded individual lines
                    ax.plot(E_F, E_f_line, color=color, lw=0.7, alpha=0.25, ls='--')

            if lines:
                lower_envelope = np.min(np.array(lines), axis=0)
                ax.plot(E_F, lower_envelope, color=color, lw=2.5, label=name)

        # Band edges
        ax.axvspan(-0.5, 0, alpha=0.12, color='cornflowerblue')
        ax.axvspan(band_gap, band_gap + 0.5, alpha=0.12, color='orange')
        ax.axhline(0, color='k', ls='--', alpha=0.3, lw=0.8)

        ax.set_xlim(-0.2, band_gap + 0.2)
        y_min = min(0, min(env.min() for env in [np.min(np.array(
            [ld[limit_name] + q * E_F for q, ld in qd.items() if limit_name in ld]),
            axis=0) for qd in all_results.values() if qd]) - 0.5)
        ax.set_ylim(y_min, None)
        ax.set_xlabel("Fermi Level (eV)", fontsize=14)
        ax.set_ylabel("Formation Energy (eV)", fontsize=14)
        ax.set_title(f"Formation Energy Diagram ({limit_name})", fontsize=14)
        ax.legend(fontsize=11, loc='best')
        plt.tight_layout()

        fig_path = os.path.join(args.output, f'formation_energy_{limit_name}.pdf')
        fig.savefig(fig_path, dpi=300)
        fig_path_png = os.path.join(args.output, f'formation_energy_{limit_name}.png')
        fig.savefig(fig_path_png, dpi=150)
        plt.close(fig)
        log(f"\nPlot saved: {fig_path}")

    # ---- 保存日志 ----
    with open(log_path, 'w') as f:
        f.write('\n'.join(log_lines) + '\n')
    log(f"Log saved: {log_path}")

if __name__ == '__main__':
    main()

formation_config.yaml 配置文件模板

# 形成能计算配置
bulk_static: Bulk/statics