功函数

#!/usr/bin/env python3
"""Rebuild planar potential and vacuum-referenced levels from VASP outputs.

Required files in the current directory: LOCPOT, OUTCAR, EIGENVAL.
The reader is intentionally restricted to a non-spin-polarized, gapped,
single-scalar LOCPOT case. It rejects unsupported inputs instead of guessing.
"""
from __future__ import annotations

import argparse
import hashlib
import json
import math
import re
from pathlib import Path

import numpy as np

A = np.asarray


def sha256(path: Path) -> str:
    h = hashlib.sha256()
    with path.open("rb") as stream:
        for block in iter(lambda: stream.read(1024 * 1024), b""):
            h.update(block)
    return h.hexdigest()


def read_locpot(path: Path):
    lines = path.read_text(errors="strict").splitlines()
    if len(lines) < 10:
        raise ValueError("LOCPOT is too short to contain a POSCAR header and grid")
    scale = float(lines[1].split()[0])
    if scale <= 0:
        raise ValueError("This example reader requires a positive POSCAR scale")
    cell = A([[float(x) * scale for x in lines[i].split()[:3]]
              for i in (2, 3, 4)], dtype=float)
    index = 5
    species_or_counts = lines[index].split()
    index += 1
    if all(re.fullmatch(r"\d+", word) for word in species_or_counts):
        counts = [int(word) for word in species_or_counts]
    else:
        counts = [int(word) for word in lines[index].split()]
        index += 1
    mode = lines[index].strip().lower()
    index += 1
    if mode.startswith("s"):
        mode = lines[index].strip().lower()
        index += 1
    if not (mode.startswith("d") or mode.startswith("c") or mode.startswith("k")):
        raise ValueError(f"Unrecognized coordinate mode in LOCPOT header: {mode!r}")
    natoms = sum(counts)
    if natoms <= 0:
        raise ValueError("Invalid atom counts in LOCPOT header")
    index += natoms
    while index < len(lines) and not lines[index].strip():
        index += 1
    grid = tuple(map(int, lines[index].split()))
    index += 1
    if len(grid) != 3 or min(grid) <= 0:
        raise ValueError(f"Invalid LOCPOT grid dimensions: {grid}")
    nx, ny, nz = grid
    expected = nx * ny * nz
    fields = []
    for line in lines[index:]:
        fields.extend(float(token.replace("D", "E").replace("d", "e"))
                      for token in line.split())
    if len(fields) != expected:
        raise ValueError(
            f"Expected one scalar potential block ({expected} values), found {len(fields)}"
        )
    values = np.asarray(fields, dtype=float)
    if not np.isfinite(values).all():
        raise ValueError("LOCPOT contains non-finite potential values")
    normal = np.cross(cell[0], cell[1])
    area = np.linalg.norm(normal)
    if area == 0:
        raise ValueError("The first two lattice vectors do not span a surface")
    normal_height = abs(float(np.dot(cell[2], normal))) / area
    if normal_height <= 0:
        raise ValueError("Invalid cell height along the surface normal")
    # VASP writes x fastest, then y, with z as the slowest index.
    field = values.reshape((nz, ny, nx))
    planar = field.mean(axis=(1, 2))
    z = np.arange(nz, dtype=float) * normal_height / nz
    return grid, normal_height, z, planar


def outcar_values(path: Path):
    text = path.read_text(errors="strict")
    if "aborting loop because EDIFF is reached" not in text:
        raise ValueError("OUTCAR has no EDIFF convergence marker")
    if "General timing and accounting" not in text:
        raise ValueError("OUTCAR has no final timing/accounting section")
    if "LVHAR" not in text or not re.search(r"LVHAR\s*=\s*T\b", text):
        raise ValueError("OUTCAR does not confirm LVHAR=T")
    if re.search(r"LVTOT\s*=\s*T\b", text):
        raise ValueError("OUTCAR also has LVTOT=T; this route expects LVHAR only")
    ispin = re.findall(r"^\s*ISPIN\s*=\s*(\d+)", text, flags=re.M)
    nelect = re.findall(r"^\s*NELECT\s*=\s*([-+0-9.]+)", text, flags=re.M)
    efermi = re.findall(r"E-fermi\s*:\s*([-+0-9.]+)", text)
    if not ispin or int(ispin[-1]) != 1:
        raise ValueError("This band-edge reader requires a non-spin-polarized ISPIN=1 run")
    if not nelect or not efermi:
        raise ValueError("OUTCAR is missing NELECT or E-fermi")
    return int(round(float(nelect[-1]))), float(efermi[-1])


def read_eigenval(path: Path, expected_electrons: int):
    with path.open() as stream:
        for _ in range(5):
            if not stream.readline():
                raise ValueError("EIGENVAL ended before its electron/k-point/band header")
        try:
            electrons, nkpoints, nbands = map(int, stream.readline().split())
        except Exception as exc:
            raise ValueError("Invalid EIGENVAL electron/k-point/band header") from exc
        if electrons != expected_electrons:
            raise ValueError(
                f"EIGENVAL has {electrons} electrons but OUTCAR has {expected_electrons}"
            )
        if electrons % 2:
            raise ValueError("An even electron count is required for this ISPIN=1 example")
        nocc = electrons // 2
        if not (0 < nocc < nbands):
            raise ValueError("EIGENVAL does not contain both occupied and empty bands")
        k_weights, energies, occupations = [], [], []
        for _ in range(nkpoints):
            line = stream.readline()
            while line and not line.strip():
                line = stream.readline()
            if not line:
                raise ValueError("EIGENVAL ended before all k-point blocks were read")
            point = list(map(float, line.split()))
            if len(point) != 4:
                raise ValueError("Expected kx ky kz weight on each EIGENVAL k-point line")
            k_weights.append(point[3])
            e_k, occ_k = [], []
            for _ in range(nbands):
                row = stream.readline().split()
                if len(row) != 3:
                    raise ValueError("Expected band index, energy, and occupation (ISPIN=1)")
                band_index, energy, occupation = map(float, row)
                e_k.append(energy)
                occ_k.append(occupation)
            energies.append(e_k)
            occupations.append(occ_k)
    weights = np.asarray(k_weights, dtype=float)
    eigenvalues = np.asarray(energies, dtype=float)
    occ = np.asarray(occupations, dtype=float)
    if abs(float(weights.sum()) - 1.0) > 1e-5:
        raise ValueError(f"EIGENVAL k-point weights sum to {weights.sum():.8g}, not 1")
    weighted_electrons = 2.0 * float(np.sum(weights[:, None] * occ))
    if abs(weighted_electrons - expected_electrons) > 1e-3:
        raise ValueError(
            f"Weighted EIGENVAL occupations give {weighted_electrons:.8f} electrons, "
            f"not {expected_electrons}"
        )
    if float(occ[:, nocc - 1].min()) < 0.5 or float(occ[:, nocc].max()) > 0.5:
        raise ValueError(
            "The NELECT/2 band boundary is partially occupied; this is not a gapped "
            "non-spin-polarized case"
        )
    vbm = float(eigenvalues[:, nocc - 1].max())
    cbm = float(eigenvalues[:, nocc].min())
    if cbm <= vbm:
        raise ValueError("Sampled EIGENVAL band edges do not form a positive gap")
    return {
        "electrons": electrons,
        "nkpoints": nkpoints,
        "bands": nbands,
        "weighted_electrons": weighted_electrons,
        "vbm_eV": vbm,
        "cbm_eV": cbm,
        "indirect_gap_eV": cbm - vbm,
    }


def parse_window(spec: str):
    try:
        low, high = map(float, spec.split(":"))
    except Exception as exc:
        raise argparse.ArgumentTypeError("Use a window such as 1:3 (angstrom)") from exc
    if not (math.isfinite(low) and math.isfinite(high) and low < high):
        raise argparse.ArgumentTypeError("Window endpoints must be finite with low < high")
    return low, high


def main():
    parser = argparse.ArgumentParser(
        description="Average LVHAR from LOCPOT and align the gapped band edges to vacuum."
    )
    parser.add_argument(
        "--windows", nargs="+", type=parse_window, default=[(1.0, 3.0), (15.0, 17.0)],
        metavar="LOW:HIGH", help="vacuum windows in angstrom (default: 1:3 15:17)"
    )
    args = parser.parse_args()
    if not args.windows:
        raise ValueError("At least one vacuum window is required")
    loct = Path("LOCPOT")
    outcar = Path("OUTCAR")
    eigenval = Path("EIGENVAL")
    if not all(path.is_file() for path in (loct, outcar, eigenval)):
        raise FileNotFoundError("Run in a directory containing LOCPOT, OUTCAR, EIGENVAL")

    grid, height, z, potential = read_locpot(loct)
    nelect, ef = outcar_values(outcar)
    bands = read_eigenval(eigenval, nelect)
    windows = []
    for low, high in args.windows:
        if low < 0 or high > height:
            raise ValueError(
                f"Window {low:g}:{high:g} A lies outside 0:{height:.8f} A"
            )
        chosen = (z >= low) & (z <= high)
        values = potential[chosen]
        if values.size == 0:
            raise ValueError(f"Window {low:g}:{high:g} A contains no LOCPOT planes")
        mean = float(values.mean())
        row = {
            "lo_A": low, "hi_A": high, "n": int(values.size),
            "mean_eV": mean,
            "std_eV": float(values.std(ddof=0)),
            "range_eV": float(values.max() - values.min()),
            "vacuum_minus_fermi_eV": mean - ef,
            "vacuum_minus_vbm_eV": mean - bands["vbm_eV"],
            "vacuum_minus_cbm_eV": mean - bands["cbm_eV"],
        }
        windows.append(row)

    order = sorted(windows, key=lambda item: item["lo_A"])
    for left, right in zip(order, order[1:]):
        if left["hi_A"] > right["lo_A"]:
            raise ValueError("Vacuum windows overlap")

    np.savetxt(
        "PLANAR_AVERAGE.dat", np.column_stack((z, potential)),
        fmt=("%.10f", "%.12f"),
        header="z_A  planar_potential_eV",
        comments="# ",
    )
    result = {
        "code": "VASP",
        "potential_component": "LVHAR (ionic + Hartree; LVTOT is false)",
        "surface_normal": "normal to lattice vectors a and b",
        "normal_height_A": height,
        "grid": list(grid),
        "scalar_values": int(np.prod(grid)),
        "ispin": 1,
        "electrons": nelect,
        "nkpoints": bands["nkpoints"],
        "bands": bands["bands"],
        "weighted_electrons": bands["weighted_electrons"],
        "fermi_eV": ef,
        "vbm_eV": bands["vbm_eV"],
        "cbm_eV": bands["cbm_eV"],
        "indirect_gap_eV": bands["indirect_gap_eV"],
        "formula": {
            "work_function_eV": "V_vacuum - E_F",
            "ionization_potential_eV": "V_vacuum - VBM",
            "electron_affinity_eV": "V_vacuum - CBM",
            "potential_spread_eV": "max(Vbar) - min(Vbar) within each selected window",
        },
        "windows": windows,
        "input_sha256": {
            path.name: sha256(path) for path in (loct, outcar, eigenval)
        },
        "limits": [
            "EIGENVAL band edges are extrema on the sampled SCF k mesh, not a continuous-Brillouin-zone search.",
            "E_F in a semiconductor is the chemical potential printed for this occupation setup; it can move within the gap.",
            "Window spread is a local flatness diagnostic, not an uncertainty or convergence estimate.",
        ],
    }
    Path("workfunction-summary.json").write_text(
        json.dumps(result, indent=2, sort_keys=True) + "\n"
    )
    print(
        f"grid={grid}; scalar values={int(np.prod(grid))}; normal height={height:.10f} A"
    )
    print(
        f"NELECT={nelect}; NKPTS={bands['nkpoints']}; weighted electrons="
        f"{bands['weighted_electrons']:.8f}; E_F={ef:.6f} eV"
    )
    print(
        f"sampled VBM={bands['vbm_eV']:.6f} eV; CBM={bands['cbm_eV']:.6f} eV; "
        f"gap={bands['indirect_gap_eV']:.6f} eV"
    )
    for item in windows:
        print(
            f"z={item['lo_A']:.2f}:{item['hi_A']:.2f} A N={item['n']} "
            f"V_vac={item['mean_eV']:.9f} eV std={item['std_eV']:.6g} eV "
            f"range={item['range_eV']:.6g} eV Phi={item['vacuum_minus_fermi_eV']:.9f} eV"
        )


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Plot the actual full-cell planar LVHAR profile used to inspect the vacuum plateau."""
from __future__ import annotations

import json
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.patches import Patch
import numpy as np

from atlas_plot_style import install

install()


def poscar_atomic_extent(path: Path, normal_height: float, normal: np.ndarray) -> tuple[float, float]:
    """Return the z extent of atomic coordinates projected along the cell normal."""
    lines = path.read_text().splitlines()
    scale = float(lines[1].split()[0])
    cell = np.asarray([[float(x) for x in lines[i].split()[:3]] for i in (2, 3, 4)]) * scale
    i = 5
    fields = lines[i].split()
    i += 1
    if all(word.lstrip("+").isdigit() for word in fields):
        counts = [int(word) for word in fields]
    else:
        counts = [int(word) for word in lines[i].split()]
        i += 1
    if lines[i].strip().lower().startswith("s"):
        i += 1
    mode = lines[i].strip().lower()
    i += 1
    coordinates = np.asarray([
        [float(x) for x in lines[i + j].split()[:3]]
        for j in range(sum(counts))
    ])
    if mode.startswith("d"):
        cartesian = coordinates @ cell
    elif mode.startswith(("c", "k")):
        cartesian = coordinates * scale
    else:
        raise ValueError(f"Unsupported POSCAR coordinate mode: {mode!r}")
    z = np.mod(cartesian @ normal, normal_height)
    return float(z.min()), float(z.max())


profile_path = Path("PLANAR_AVERAGE.dat")
summary_path = Path("workfunction-summary.json")
if not profile_path.is_file() or not summary_path.is_file():
    raise FileNotFoundError("Run analyze_workfunction.py first")

profile = np.loadtxt(profile_path, comments="#")
summary = json.loads(summary_path.read_text())
windows = sorted(summary["windows"], key=lambda item: item["lo_A"])
if profile.ndim != 2 or profile.shape[1] != 2:
    raise ValueError("Expected z_A and planar_potential_eV columns")
if len(windows) != 2:
    raise ValueError("This example expects its two measured surface vacuum windows")
z, potential = profile[:, 0], profile[:, 1]
height = float(summary["normal_height_A"])
if not np.all(np.isfinite(profile)) or z.max() >= height:
    raise ValueError("Profile contains invalid values or an out-of-cell z coordinate")

cell = np.asarray([[float(x) for x in line.split()[:3]]
                   for line in Path("POSCAR").read_text().splitlines()[2:5]], dtype=float)
scale = float(Path("POSCAR").read_text().splitlines()[1].split()[0])
cell *= scale
normal = np.cross(cell[0], cell[1])
normal /= np.linalg.norm(normal)
if np.dot(cell[2], normal) < 0:
    normal *= -1.0
atom_lo, atom_hi = poscar_atomic_extent(Path("POSCAR"), height, normal)

vacuum_reference = float(windows[0]["mean_eV"])
relative_potential = potential - vacuum_reference
fermi_relative = float(summary["fermi_eV"]) - vacuum_reference
work_function = vacuum_reference - float(summary["fermi_eV"])
colors = ("#0072b2", "#d55e00")

fig, ax = plt.subplots(figsize=(8.1, 4.8), layout="constrained")
ax.plot(z, relative_potential, color="#222222", lw=1.25,
        label=r"Planar-averaged $V_{\mathrm{LVHAR}}(z)$")
ax.axhline(0.0, color=colors[0], ls="--", lw=1.0,
           label="Lower-z vacuum reference")
ax.axhline(fermi_relative, color=colors[1], ls="-.", lw=1.0,
           label=r"$E_F$ from the same SCF")
ax.axvspan(atom_lo, atom_hi, color="#777777", alpha=0.12, zorder=0)
for window in windows:
    ax.axvspan(window["lo_A"], window["hi_A"], color=colors[0], alpha=0.08, zorder=0)

handles = [
    Line2D([0], [0], color="#222222", lw=1.25,
           label=r"Planar-averaged $V_{\mathrm{LVHAR}}(z)$"),
    Line2D([0], [0], color=colors[0], ls="--", lw=1.0,
           label="Lower-z vacuum reference"),
    Line2D([0], [0], color=colors[1], ls="-.", lw=1.0,
           label=r"$E_F$ from the same SCF"),
    Patch(facecolor="#777777", alpha=0.12, label="Atomic z extent from POSCAR"),
    Patch(facecolor=colors[0], alpha=0.08, label="Vacuum windows used in the table"),
]
ax.legend(handles=handles, frameon=False, ncols=2, loc="lower left")
arrow_x = 2.55
ax.annotate(
    "", xy=(arrow_x, fermi_relative), xytext=(arrow_x, 0.0),
    arrowprops={"arrowstyle": "<->", "color": colors[1], "lw": 1.1},
)
ax.text(arrow_x + 0.22, 0.5 * fermi_relative,
        rf"$\Phi={work_function:.4f}\ \mathrm{{eV}}$",
        color=colors[1], ha="left", va="center")

span = float(relative_potential.max() - relative_potential.min())
ax.set_xlim(0.0, height)
ax.set_ylim(relative_potential.min() - 0.04 * span,
            relative_potential.max() + 0.04 * span)
ax.set_xlabel("Distance along the surface normal (Å)")
ax.set_ylabel(r"$\bar{V}(z)-V_{\mathrm{vac,lower}}$ (eV)")
ax.set_title(r"SnSe$_2$ monolayer: planar-averaged LVHAR and work-function reference")
ax.grid(False)

output = "interface-magnet-workfunction-profile.png"
fig.savefig(output)
plt.close(fig)

upper_minus_lower = float(windows[1]["mean_eV"] - windows[0]["mean_eV"])
print(f"Wrote interface-magnet-workfunction-profile.png, .svg, and .pdf")
print(f"Lower-z Vvac = {vacuum_reference:.9f} eV; E_F = {summary['fermi_eV']:.6f} eV")
print(f"Phi at the printed E_F = {work_function:.9f} eV")
print(f"Upper-minus-lower vacuum mean = {upper_minus_lower * 1000.0:+.5f} meV")
print("That last difference is reported numerically only; it is not resolved as a material effect.")

from __future__ import print_function
import sys, math, json, hashlib

def cross(a,b):
    return [a[1]*b[2]-a[2]*b[1],a[2]*b[0]-a[0]*b[2],a[0]*b[1]-a[1]*b[0]]
def dot(a,b): return sum(x*y for x,y in zip(a,b))

def read_grid(name):
    f=open(name)
    title=f.readline().strip()
    scale=float(f.readline().split()[0])
    cell=[[float(x)*scale for x in f.readline().split()[:3]] for i in range(3)]
    if scale <= 0: raise ValueError('This reader requires a positive POSCAR scale')
    words=f.readline().split()
    if all(x.isdigit() for x in words):
        counts=list(map(int,words))
    else:
        counts=list(map(int,f.readline().split()))
    mode=f.readline().strip()
    if mode.lower().startswith('s'): mode=f.readline().strip()
    coords=[f.readline().split()[:3] for i in range(sum(counts))]
    line=f.readline()
    while line and not line.strip(): line=f.readline()
    grid=list(map(int,line.split()))
    if len(grid)!=3 or min(grid)<=0: raise ValueError('Invalid FFT grid')
    n=grid[0]*grid[1]*grid[2]
    vals=[]
    while len(vals)<n:
        line=f.readline()
        if not line: raise ValueError('Truncated potential: %d/%d'%(len(vals),n))
        vals.extend(float(x.replace('D','E')) for x in line.split())
    if len(vals)!=n: raise ValueError('Unexpected extra values in scalar block')
    f.close()
    if any(math.isnan(x) or math.isinf(x) for x in vals):
        raise ValueError("Non-finite potential value")
    return cell,grid,vals

if __name__=='__main__':
    name=sys.argv[1] if len(sys.argv)>1 else 'LOCPOT'
    cell,grid,v=read_grid(name)
    area=math.sqrt(dot(cross(cell[0],cell[1]),cross(cell[0],cell[1])))
    height=abs(dot(cell[2],cross(cell[0],cell[1])))/area
    nxy=grid[0]*grid[1]
    avg=[sum(v[i*nxy:(i+1)*nxy])/nxy for i in range(grid[2])]
    zz=[height*i/grid[2] for i in range(grid[2])]
    with open('PLANAR_AVERAGE.dat','w') as f:
        f.write('# z_A  planar_potential_eV\n')
        for z,p in zip(zz,avg): f.write('%.10f %.12f\n'%(z,p))
    summary={'grid':grid,'points':len(v),'normal_height_A':height,'source_sha256':hashlib.sha256(open(name,'rb').read()).hexdigest(),'windows':[]}
    print('grid = %d %d %d; scalar values = %d'%tuple(grid+[len(v)]))
    print('normal height = %.10f A; output = PLANAR_AVERAGE.dat'%height)
    for spec in sys.argv[2:]:
        lo,hi=map(float,spec.split(':'))
        a=[p for z,p in zip(zz,avg) if lo<=z<=hi]
        if not a: raise ValueError('Empty averaging window')
        mean=sum(a)/len(a); std=math.sqrt(sum((x-mean)**2 for x in a)/len(a)); span=max(a)-min(a)
        row={'lo_A':lo,'hi_A':hi,'n':len(a),'mean_eV':mean,'std_eV':std,'range_eV':span}
        summary['windows'].append(row)
        print('window %.2f:%.2f A  N=%d  mean=%.9f eV  std=%.6g eV  range=%.6g eV'%(lo,hi,len(a),mean,std,span))
    with open('potential-summary.json','w') as f: json.dump(summary,f,indent=2,sort_keys=True)

from __future__ import print_function
import json,re,hashlib
out=open('OUTCAR').read()
if 'aborting loop because EDIFF is reached' not in out:
    raise ValueError('Electronic convergence line is absent')
if 'General timing and accounting' not in out:
    raise ValueError('Normal final accounting section is absent')
ef=float(re.findall(r'E-fermi\s*:\s*([-+0-9.]+)',out)[-1])
p=json.load(open('potential-summary.json'))
with open('EIGENVAL') as f:
    for i in range(5): f.readline()
    ne,nk,nb=map(int,f.readline().split())
    if ne % 2: raise ValueError('This band-edge reader expects even-electron, non-spin-polarized input')
    occupied=[];empty=[]
    for ik in range(nk):
        line=f.readline()
        while line and not line.strip(): line=f.readline()
        if not line: raise ValueError('Truncated EIGENVAL before k point')
        if len(line.split())!=4: raise ValueError('Invalid k-point line')
        bands=[list(map(float,f.readline().split())) for ib in range(nb)]
        if any(len(row)!=3 for row in bands): raise ValueError('Invalid or truncated non-spin EIGENVAL band block')
        occupied.append(bands[ne//2-1][1]);empty.append(bands[ne//2][1])
vbm=max(occupied);cbm=min(empty)
r={'fermi_eV':ef,'vbm_eV':vbm,'cbm_eV':cbm,'indirect_gap_eV':cbm-vbm,'nkpoints':nk,'bands':nb,'electrons':ne,'windows':[]}
for w in p['windows']:
    item=dict(w)
    item['vacuum_minus_fermi_eV']=w['mean_eV']-ef
    item['vacuum_minus_vbm_eV']=w['mean_eV']-vbm
    item['vacuum_minus_cbm_eV']=w['mean_eV']-cbm
    r['windows'].append(item)
json.dump(r,open('workfunction-summary.json','w'),indent=2)
print('E_F = %.6f eV; VBM = %.6f eV; CBM = %.6f eV; gap = %.6f eV'%(ef,vbm,cbm,cbm-vbm))
for w in r['windows']:
    print('z = %.2f:%.2f A; V_vac = %.9f eV; Phi(E_F) = %.9f eV; V_vac-VBM = %.9f eV; V_vac-CBM = %.9f eV'%(w['lo_A'],w['hi_A'],w['mean_eV'],w['vacuum_minus_fermi_eV'],w['vacuum_minus_vbm_eV'],w['vacuum_minus_cbm_eV']))

[bcgong@localhost snse2_workfunction]$ tail -8 OSZICAR
DAV:  16    -0.118616121789E+02   -0.11011E-04   -0.16797E-07  4460   0.284E-03    0.893E-04
DAV:  17    -0.118616188163E+02   -0.66374E-05   -0.38095E-07  4636   0.216E-03    0.126E-03
DAV:  18    -0.118616239345E+02   -0.51182E-05   -0.15024E-07  4476   0.217E-03    0.320E-04
DAV:  19    -0.118616257043E+02   -0.17698E-05   -0.56160E-08  4528   0.662E-04    0.339E-04
DAV:  20    -0.118616260960E+02   -0.39175E-06   -0.97347E-09  4376   0.391E-04    0.116E-04
DAV:  21    -0.118616262047E+02   -0.10867E-06   -0.18647E-09  4340   0.150E-04    0.186E-05
DAV:  22    -0.118616262492E+02   -0.44434E-07   -0.62584E-10  3884   0.883E-05
   1 F= -.12268138E+02 E0= -.12268138E+02  d E =-.349964E-11

[bcgong@localhost snse2_workfunction]$ head -7 OUTCAR
 vasp.5.4.4.18Apr17-6-g9f103f2a35 (build Feb 26 2024 21:30:50) complex          
  
 executed on             LinuxIFC date 2026.09.22  22:13:18
 running on    8 total cores
 distrk:  each k-point on    8 cores,    1 groups
 distr:  one band on NCORES_PER_BAND=   2 cores,    4 groups

[bcgong@localhost snse2_workfunction]$ ls -lh OUTCAR OSZICAR CHGCAR LOCPOT EIGENVAL WAVECAR
-rw-rw-r-- 1 bcgong bcgong  18M Sep 22 22:14 CHGCAR
-rw-rw-r-- 1 bcgong bcgong  77K Sep 22 22:14 EIGENVAL
-rw-rw-r-- 1 bcgong bcgong  18M Sep 22 22:14 LOCPOT
-rw-rw-r-- 1 bcgong bcgong 2.1K Sep 22 22:14 OSZICAR
-rw-rw-r-- 1 bcgong bcgong 189K Sep 22 22:14 OUTCAR
-rw-rw-r-- 1 bcgong bcgong    0 Sep 22 22:13 WAVECAR

[bcgong@localhost snse2_workfunction]$ head -5 PLANAR_AVERAGE.dat
# z_A  planar_potential_eV
0.0000000000 3.306285072550
0.0655617786 3.306304853416
0.1311235571 3.306285038010
0.1966853357 3.306304820274

固定结构 + 同一套赝势 / 均匀 k 网格
  └─ SCF + LVHAR
       ├─ OUTCAR → 收敛与 E_F
       ├─ EIGENVAL → 当前网格的 VBM / CBM
       └─ LOCPOT → 平面平均 → 平坦真空窗口
                                    └─ V_vac − E_F；同时注明带隙中的化学势