占据态 Berry 相位与切片陈数

"""FHS determinant links from native QE/Wannier90 overlaps; NumPy only."""
from pathlib import Path
import csv
import hashlib
import json
import re
import time
import xml.etree.ElementTree as ET
import numpy as np

ROOT = Path(__file__).resolve().parent
OUT = ROOT / "results"
OUT.mkdir(exist_ok=True)

def block(text, name):
    return re.search(r"begin\s+"+name+r"\s*\n(.*?)end\s+"+name,
                     text, re.S | re.I).group(1).strip().splitlines()

def write_csv(name, rows):
    with (OUT / name).open("w") as stream:
        writer = csv.DictWriter(stream, fieldnames=rows[0].keys())
        writer.writeheader()
        writer.writerows(rows)

def read_case(n):
    directory = ROOT / "source" / f"k{n}"
    win = (directory / "silicon.win").read_text()
    points = np.array([list(map(float, line.split())) for line in block(win, "kpoints")])
    indices = np.rint(points * n).astype(int)
    assert np.max(np.abs(points * n - indices)) < 1e-8
    assert len(set(map(tuple, indices))) == n**3
    lookup = {tuple(v): i for i, v in enumerate(indices)}
    nnkp = (directory / "silicon.nnkp").read_text()
    nn_points = block(nnkp, "kpoints")
    assert int(nn_points[0]) == len(points)
    assert np.allclose(points, np.array([list(map(float, x.split())) for x in nn_points[1:]]), atol=1e-8, rtol=0)
    bvecs = np.array([list(map(float, line.split())) for line in block(nnkp, "recip_lattice")])
    real = np.array([list(map(float, line.split())) for line in block(nnkp, "real_lattice")])
    reciprocal_error = float(np.max(np.abs(real @ bvecs.T - 2*np.pi*np.eye(3))))
    assert reciprocal_error < 1e-6
    nn_lines = block(nnkp, "nnkpts")
    nnb = int(nn_lines[0])
    declared = {tuple(map(int, line.split())) for line in nn_lines[1:]}
    with (directory / "silicon.mmn").open() as stream:
        stream.readline()
        nb, nk, nn = map(int, stream.readline().split())
        assert (nb, nk, nn) == (4, n**3, nnb)
        matrices = {}
        for _ in range(nk*nn):
            header = tuple(map(int, stream.readline().split()))
            assert len(header) == 5 and header not in matrices
            values = [complex(*map(float, stream.readline().split())) for _ in range(nb*nb)]
            matrices[header] = np.array(values).reshape((nb, nb), order="F")
        assert not stream.read().strip()
    assert set(matrices) == declared
    tree = ET.parse(directory / "nscf.data-file-schema.xml").getroot()
    vals = lambda tag: [e.text.strip() for e in tree.iter() if e.tag.split('}')[-1] == tag]
    assert set(vals("nbnd")) == {"4"} and set(vals("nks")) == {str(n**3)}
    assert all(float(x) == 8 for x in vals("nelec"))
    for tag in ["lsda", "noncolin", "spinorbit"]:
        assert set(vals(tag)) == {"false"}
    occ = [e for e in vals("occupations") if e != "fixed"]
    assert len(occ) == n**3 and all(np.array_equal(np.fromstring(x, sep=" "), np.ones(4)) for x in occ)
    for stem in ["si.scf", "si.nscf", "pw2wan"]:
        text = (directory / (stem+".out")).read_text()
        assert text.count("JOB DONE.") == 1
        assert not re.search(r"Error in routine|convergence NOT|eigenvalues not converged|MPI_ABORT", text)
        assert (directory / (stem+".err")).stat().st_size == 0
    selected = {}
    reverse_error = 0.0
    for header, matrix in matrices.items():
        i, j, *g = header
        reverse = (j, i, *[-v for v in g])
        assert reverse in matrices
        reverse_error = max(reverse_error, float(np.max(np.abs(matrix-matrices[reverse].conj().T))))
        delta = (points[j-1]+g-points[i-1])*n
        step = np.rint(delta).astype(int)
        assert np.max(np.abs(delta-step)) < 1e-8
        for axis in [0, 1]:
            unit = np.eye(3, dtype=int)[axis]
            if np.array_equal(step, unit):
                assert (i-1, axis) not in selected
                assert j-1 == lookup[tuple((indices[i-1]+unit) % n)]
                selected[i-1, axis] = (j-1, tuple(g), matrix)
    assert len(selected) == 2*n**3 and reverse_error < 1e-9
    return directory, points, indices, lookup, selected, bvecs, reciprocal_error, reverse_error

def links_and_flux(n, indices, lookup, selected, gauge=None):
    links = {}
    for (i, axis), (j, g, matrix) in selected.items():
        if gauge is not None:
            matrix = gauge[i].conj().T @ matrix @ gauge[j]
        determinant = np.linalg.det(matrix)
        assert abs(determinant) > 1e-12
        links[i, axis] = determinant / abs(determinant)
    flux = np.empty(n**3)
    for i, coord in enumerate(indices):
        i1 = lookup[tuple((coord + [1, 0, 0]) % n)]
        i2 = lookup[tuple((coord + [0, 1, 0]) % n)]
        flux[i] = np.angle(links[i,0]*links[i1,1]*np.conj(links[i2,0])*np.conj(links[i,1]))
    return links, flux

def run(n):
    start = time.perf_counter()
    directory, points, indices, lookup, selected, bvecs, reciprocal_error, reverse_error = read_case(n)
    area = float(np.linalg.norm(np.cross(bvecs[0], bvecs[1])) / n**2)
    rows = []
    for (i, axis), (j, g, matrix) in selected.items():
        sv = np.linalg.svd(matrix, compute_uv=False)
        rows.append(dict(grid=n,k_index=i+1,axis=axis+1,neighbor_index=j+1,G1=g[0],G2=g[1],G3=g[2],min_singular=float(sv[-1]),max_singular=float(sv[0]),abs_determinant=float(abs(np.linalg.det(matrix)))))
    minsv = min(x["min_singular"] for x in rows)
    maxsv = max(x["max_singular"] for x in rows)
    assert minsv > 1e-8 and maxsv < 1.001
    links, flux = links_and_flux(n, indices, lookup, selected)
    link_norm_error = max(abs(abs(v)-1) for v in links.values())
    assert link_norm_error < 1e-14
    gauge_rows = []
    for trial in range(12):
        seed = 2026092200 + 100*n + trial
        rng = np.random.default_rng(seed)
        gauge = []
        for _ in range(n**3):
            z = rng.normal(size=(4,4)) + 1j*rng.normal(size=(4,4))
            q, rr = np.linalg.qr(z)
            q = q @ np.diag(np.diag(rr)/np.abs(np.diag(rr)))
            gauge.append(q)
        gauge = np.array(gauge)
        norm_error = float(np.max(np.abs(gauge.conj().transpose(0,2,1)@gauge-np.eye(4))))
        _, changed = links_and_flux(n, indices, lookup, selected, gauge)
        flux_error = float(np.max(np.abs(np.angle(np.exp(1j*(changed-flux))))))
        chern_error = max(abs(np.sum(changed[indices[:,2]==s]-flux[indices[:,2]==s])/(2*np.pi)) for s in range(n))
        assert norm_error < 1e-12 and flux_error < 1e-12 and chern_error < 1e-12
        gauge_rows.append(dict(grid=n,trial=trial,seed=seed,unitary_error=norm_error,max_flux_difference_rad=flux_error,max_chern_difference=float(chern_error)))
    slice_rows=[]
    for s in range(n):
        f = flux[indices[:,2]==s]
        chern = float(np.sum(f)/(2*np.pi))
        integer = int(np.rint(chern))
        assert abs(chern-integer) < 1e-10
        slice_rows.append(dict(grid=n,k3_index=s,k3_fraction=s/n,plaquettes=n*n,chern_raw=chern,chern_integer=integer,max_abs_phase_rad=float(np.max(np.abs(f))),plaquette_area_invA2=area))
    write_csv(f"k{n}-links.csv", rows)
    write_csv(f"k{n}-gauge-check.csv", gauge_rows)
    write_csv(f"k{n}-plaquettes.csv", [dict(grid=n,k_index=i+1,k1_index=int(c[0]),k2_index=int(c[1]),k3_index=int(c[2]),k1_fraction=float(points[i,0]),k2_fraction=float(points[i,1]),k3_fraction=float(points[i,2]),phase_rad=float(flux[i]),phase_per_area_A2=float(flux[i]/area)) for i,c in enumerate(indices)])
    summary=dict(grid=n,occupied_bands=4,kpoints=n**3,mmn_neighbors=8,selected_links=len(selected),boundary_links=sum(any(x['G'+str(i)] for i in [1,2,3]) for x in rows),min_link_singular_value=minsv,max_link_singular_value=maxsv,min_abs_determinant=min(x['abs_determinant'] for x in rows),reciprocal_duality_error=reciprocal_error,reverse_link_max_error=reverse_error,normalized_link_norm_error=link_norm_error,max_abs_phase_rad=float(np.max(np.abs(flux))),branch_margin_rad=float(np.pi-np.max(np.abs(flux))),random_gauge_trials=12,max_gauge_flux_difference_rad=max(x['max_flux_difference_rad'] for x in gauge_rows),max_gauge_chern_difference=max(x['max_chern_difference'] for x in gauge_rows),slices=slice_rows,wall_seconds=time.perf_counter()-start)
    print(f"GRID {n}x{n}x{n}: {n**3} k points, 4 occupied bands, {len(selected)} directed links")
    print(f"  singular values: min={minsv:.12f} max={maxsv:.12f}; min|det M|={summary['min_abs_determinant']:.12f}")
    print(f"  reverse overlap residual={reverse_error:.3e}; max plaquette phase={summary['max_abs_phase_rad']:.6e} rad")
    print(f"  12 random U(4) gauges: max phase difference={summary['max_gauge_flux_difference_rad']:.3e} rad")
    for row in slice_rows:
        print(f"  k3={row['k3_fraction']:.9f}: C={row['chern_raw']:+.12e}, integer={row['chern_integer']}, max|phase|={row['max_abs_phase_rad']:.6e}")
    return summary, slice_rows

summaries=[]
slices=[]
for n in [4,6]:
    result, rows = run(n)
    summaries.append(result)
    slices.extend(rows)
write_csv("slices.csv", slices)
record=dict(method="Fukui-Hatsugai-Suzuki determinant links of four occupied native QE/pw2wannier90 bands",phase_convention="Arg(U1(k) U2(k+e1) conj(U1(k+e2)) conj(U2(k)))",surface="fixed fractional k3, oriented reciprocal b1,b2 torus",spin="nonmagnetic scalar calculation: one equivalent spin channel; four spatial bands, eight electrons",claim="Discrete slice Chern calculation and internal numerical checks only; no full-zone insulating-gap or Z2 validation",numpy_version=np.__version__,cases=summaries)
(OUT/"summary.json").write_text(json.dumps(record,indent=2)+"\n")
hashes={str(f.relative_to(ROOT)):hashlib.sha256(f.read_bytes()).hexdigest() for f in sorted((ROOT/'source').rglob('*')) if f.is_file()}
(OUT/'source-sha256.json').write_text(json.dumps(hashes,indent=2)+"\n")
print("POSTPROCESS_CHECKS_PASSED; full-zone gap and material topological classification not established.")

"""Independent polar-matrix loop check, plus a nonzero synthetic link test."""
from pathlib import Path
import csv, hashlib, json
import numpy as np
root=Path(__file__).resolve().parent
hashes=json.loads((root/'results/source-sha256.json').read_text())
assert all(hashlib.sha256((root/n).read_bytes()).hexdigest()==h for n,h in hashes.items())
receipts=[]
for n in [4,6]:
    raw=(root/f'source/k{n}/silicon.mmn').read_text().splitlines()
    nb,nk,nn=map(int,raw[1].split())
    matrices={}
    for at in range(2,len(raw),nb*nb+1):
        header=tuple(map(int,raw[at].split()))
        z=np.array([complex(*map(float,line.split())) for line in raw[at+1:at+1+nb*nb]]).reshape(nb,nb,order='F')
        u,sv,vh=np.linalg.svd(z)
        matrices[header]=u@vh
    rows=list(csv.DictReader((root/f'results/k{n}-links.csv').open()))
    links={(int(x['k_index'])-1,int(x['axis'])-1):(int(x['neighbor_index'])-1,matrices[(int(x['k_index']),int(x['neighbor_index']),int(x['G1']),int(x['G2']),int(x['G3']))]) for x in rows}
    phases={}
    for k in range(nk):
        k1,m1=links[k,0];k2,m2=links[k,1]
        k12,m12=links[k1,1];k21,m21=links[k2,0]
        assert k12==k21
        phases[k]=float(np.angle(np.linalg.det(m1@m12@m21.conj().T@m2.conj().T)))
    saved=list(csv.DictReader((root/f'results/k{n}-plaquettes.csv').open()))
    error=max(abs(np.angle(np.exp(1j*(phases[int(x['k_index'])-1]-float(x['phase_rad']))))) for x in saved)
    assert error<1e-12
    receipts.append(dict(grid=n,independent_polar_loop_max_phase_error_rad=error,raw_mmn_sha256=hashlib.sha256((root/f'source/k{n}/silicon.mmn').read_bytes()).hexdigest()))
    print(f'k{n}: independent SVD polar-matrix loop phase difference = {error:.3e} rad')
# Periodic link field with a known +1 total flux; not a material calculation.
n=7
u1=np.array([[np.exp(-2j*np.pi*j/(n*n)) for j in range(n)] for i in range(n)])
u2=np.ones((n,n),complex)
for i in range(n):u2[i,-1]=np.exp(2j*np.pi*i/n)
phase=np.empty((n,n))
for i in range(n):
    for j in range(n):phase[i,j]=np.angle(u1[i,j]*u2[(i+1)%n,j]*np.conj(u1[i,(j+1)%n])*np.conj(u2[i,j]))
synthetic=float(phase.sum()/(2*np.pi))
assert abs(synthetic-1)<1e-12
print(f'Synthetic periodic link field: C = {synthetic:.12f} (expected +1; algebra check only)')
assert all(hashlib.sha256((root/n).read_bytes()).hexdigest()==h for n,h in hashes.items())
receipt=dict(source_files_unchanged=len(hashes),polar_loop_checks=receipts,synthetic_link_chern=synthetic,synthetic_scope='Algorithm orientation and periodic-boundary check, not a Si/QE result')
(root/'results/independent-check.json').write_text(json.dumps(receipt,indent=2)+'\n')
print('INDEPENDENT_CHECKS_PASSED')

"""Matrix Wilson loops of the archived scalar Si occupied subspace; no Z2 inference."""
from pathlib import Path
import argparse, csv, itertools, json, re
import xml.etree.ElementTree as ET
import numpy as np

def block(text, name):
    found = re.search(r"begin\s+" + name + r"\s*\n(.*?)end\s+" + name, text, re.S | re.I)
    if found is None:
        raise ValueError("Missing block: " + name)
    return found.group(1).strip().splitlines()

def require(condition, message):
    if not condition:
        raise ValueError(message)

def circle_error(a, b):
    return min(max(abs(((a[i]-b[j]+0.5) % 1)-0.5) for i,j in enumerate(p))
               for p in itertools.permutations(range(len(a))))

def load(data, n):
    folder = data / "source" / f"k{n}"
    points = np.array([list(map(float, x.split())) for x in block((folder/"silicon.win").read_text(), "kpoints")])
    index = np.rint(points*n).astype(int)
    require(points.shape == (n**3, 3) and np.max(abs(points*n-index)) < 1e-8, "Grid mismatch")
    lookup = {tuple(k):i for i,k in enumerate(index)}
    require(len(lookup) == n**3, "Duplicate grid points")
    nnkp = (folder/"silicon.nnkp").read_text()
    nnpoints = block(nnkp, "kpoints")
    require(int(nnpoints[0]) == n**3 and np.allclose(points, np.array([list(map(float,x.split())) for x in nnpoints[1:]]), atol=1e-8, rtol=0), "NNKP order mismatch")
    declared = {tuple(map(int,x.split())) for x in block(nnkp, "nnkpts")[1:]}
    xml = ET.parse(folder/"nscf.data-file-schema.xml").getroot()
    vals = lambda tag: [e.text.strip() for e in xml.iter() if e.tag.split("}")[-1] == tag]
    require(set(vals("nbnd")) == {"4"} and set(vals("nks")) == {str(n**3)}, "XML dimension mismatch")
    require(all(float(x)==8 for x in vals("nelec")), "Not the eight-electron Si case")
    for tag in ["lsda", "noncolin", "spinorbit"]:
        require(set(vals(tag)) == {"false"}, "Not the archived scalar case")
    occupied = [x for x in vals("occupations") if x != "fixed"]
    require(len(occupied)==n**3 and all(np.array_equal(np.fromstring(x, sep=" "), np.ones(4)) for x in occupied), "Occupation mismatch")
    matrices = {}
    with (folder/"silicon.mmn").open() as f:
        f.readline()
        nb,nk,nn = map(int,f.readline().split())
        require((nb,nk,nn)==(4,n**3,8), "MMN dimension mismatch")
        for _ in range(nk*nn):
            h=tuple(map(int,f.readline().split()))
            require(len(h)==5 and h not in matrices, "Duplicate or invalid MMN header")
            matrices[h]=np.array([complex(*map(float,f.readline().split())) for _ in range(nb*nb)]).reshape(nb,nb,order="F")
        require(not f.read().strip(), "Extra MMN records")
    require(set(matrices)==declared, "MMN/NNKP headers mismatch")
    links={}; minsv=1.; reverse=0.
    for h,m in matrices.items():
        i,j,*g=h; step=(points[j-1]+g-points[i-1])*n
        if np.allclose(step,[1,0,0],atol=1e-8,rtol=0):
            require(i-1 not in links and j-1==lookup[tuple((index[i-1]+[1,0,0]) % n)], "Bad periodic link")
            back=(j,i,*[-x for x in g]); require(back in matrices, "Missing reverse link")
            reverse=max(reverse,float(np.max(abs(m-matrices[back].conj().T))))
            u,s,vh=np.linalg.svd(m); minsv=min(minsv,float(s[-1]))
            links[i-1]=(j-1,u@vh)
    require(len(links)==n**3 and minsv>1e-8 and reverse<1e-9, "Singular or inconsistent links")
    return index, lookup, links, minsv, reverse

def run(data,n):
    index,lookup,links,minsv,reverse=load(data,n)
    rng=np.random.default_rng(2026100200+n)
    gauges=[]
    for _ in range(n**3):
        q,r=np.linalg.qr(rng.normal(size=(4,4))+1j*rng.normal(size=(4,4)))
        gauges.append(q)
    rows=[]; maxunit=0.; maxdet=0.; maxgauge=0.
    for k3 in range(n):
        for k2 in range(n):
            w=np.eye(4,dtype=complex); rotated=np.eye(4,dtype=complex); detprod=1.+0j
            for k1 in range(n):
                i=lookup[(k1,k2,k3)]; j,q=links[i]
                w=w@q; rotated=rotated@(gauges[i].conj().T@q@gauges[j]); detprod*=np.linalg.det(q)
            centres=np.sort((-np.angle(np.linalg.eigvals(w))/(2*np.pi)) % 1)
            changed=np.sort((-np.angle(np.linalg.eigvals(rotated))/(2*np.pi)) % 1)
            unit=float(np.max(abs(w.conj().T@w-np.eye(4))))
            deterror=float(abs(np.angle(np.linalg.det(w)/detprod)))
            gauge=circle_error(centres,changed)
            maxunit=max(maxunit,unit); maxdet=max(maxdet,deterror); maxgauge=max(maxgauge,gauge)
            rows.append(dict(grid=n,k2_fraction=k2/n,k3_fraction=k3/n,**{f"wcc{i+1}_mod1":float(x) for i,x in enumerate(centres)},sum_wcc_mod1=float(sum(centres)%1),unitarity_error=unit,det_phase_error_rad=deterror,gauge_spectrum_error_mod1=gauge))
    require(maxunit<1e-12 and maxdet<1e-12 and maxgauge<1e-12, "Wilson loop check failed")
    summary=dict(grid=n,loop_direction="+b1 including periodic G",loops=n*n,loop_points=n,spatial_bands=4,min_link_singular_value=minsv,reverse_overlap_error=reverse,max_unitarity_error=maxunit,max_determinant_phase_error_rad=maxdet,max_random_gauge_spectrum_error_mod1=maxgauge,random_seed=2026100200+n)
    print(f"{n}^3: {n*n} closed matrix loops; {n} links/loop; four spatial-band phases")
    print(f"  unitarity={maxunit:.3e}; determinant phase={maxdet:.3e} rad; gauge spectrum={maxgauge:.3e}")
    return rows,summary

if __name__ == "__main__":
    parser=argparse.ArgumentParser(); parser.add_argument("--data",type=Path,required=True); parser.add_argument("--output",type=Path,required=True)
    args=parser.parse_args(); args.output.mkdir(parents=True,exist_ok=True)
    rows=[]; summaries=[]
    for n in [4,6]:
        r,s=run(args.data,n); rows.extend(r); summaries.append(s)
    with (args.output/"si-wilson-loops.csv").open("w") as f:
        writer=csv.DictWriter(f,fieldnames=rows[0].keys());writer.writeheader();writer.writerows(rows)
    record=dict(quantity="-Arg(eigenvalue of product of polar MMN links)/(2*pi), modulo one",subspace="Four occupied spatial bands of eight-electron nonmagnetic scalar Si, no SOC",connection="Closed +b1 loops at each sampled fractional k2,k3",interpretation="Matrix phase extraction only; no spinful Kramers partner tracking, Z2, edge spectrum or gap validation",numpy_version=np.__version__,cases=summaries)
    (args.output/"si-wilson-summary.json").write_text(json.dumps(record,indent=2)+"\n")
    print("WILSON_MATRIX_CHECKS_PASSED; no Z2 or material topology assigned")

输入 Nk1=Nk2=81 如下；另一组只把这两个参数改成 41，两个实际输入都在下载包中。WCC 实际回路积分点数由内部自适应算法确定。

&TB_FILE
 Hrfile = 'BHZ_hr.dat'
/
&CONTROL
 BulkBand_calc = T
 WannierCenter_calc = T
 SlabSS_calc = T
/
&SYSTEM
 SOC = 1
 NumOccupied = 2
 E_FERMI = 0.0
/
&PARAMETERS
 Np = 2
 Nk1 = 81
 Nk2 = 81
 Nk3 = 1
 OmegaMin = -2.5
 OmegaMax = 2.5
 OmegaNum = 401
/
LATTICE
Angstrom
3 0 0
0 3 0
0 0 10
ATOM_POSITIONS
1
Direct
C 0 0 0
PROJECTORS
2
C s pz
SURFACE
0 0 1
1 0 0
0 1 0
KPATH_SLAB
1
-X 0 -0.5 X 0 0.5
KPATH_BULK
2
X 0.5 0 0 G 0 0 0
G 0 0 0 Y 0 0.5 0
KPLANE_BULK
0 0 0
1 0 0
0 0.5 0

官方生成器原文：

#!/bin/python3
import numpy as np
import cmath

# The Hamiltonian is 
#     ( M-Bk^2    Delta_0+A*k+  ) 
#     ( Delta_0+A*k_    -M+Bk^2 )
# where k^2=kx^2+ky^2

# Case I, QSHE with band inversion and no trivial hybridization
# Delta_0=0, M*B>0, |A|>0

# Case I, QSHE with band inversion and with trivial and nontrivial hybridization
# Delta_0=0.5, M*B>0, |A|>0


# from the kp to TB we use sustitution
# k->sin(k)
# k^2->2(1-cos(k))

# Constants
dp = np.float64
pi = np.arctan(1) * 4
zi = 1j

# Lattice constants
M = 2.0
B = 1.0
A = 1.0
Delta_0=  0.0
             

# Number of Wannier functions and R points
num_wann = 4
nrpts = 7

# R coordinates
Irvec = np.zeros((3, nrpts), dtype=int)

# Hamiltonian m,n are band indexes
HmnR = np.zeros((num_wann, num_wann, nrpts), dtype=complex)

# No of degeneracy of R point
ndegen = np.ones(nrpts, dtype=int)

# Initialization of matrices
Irvec[:, :] = 0
HmnR[:, :, :] = 0.0

# 0 0 0
ir = 0
Irvec[:, ir] = [0, 0, 0]
HmnR[0, 0, ir] = M - 4 * B
HmnR[1, 1, ir] = -M + 4 * B
HmnR[2, 2, ir] = M - 4 * B
HmnR[3, 3, ir] = -M + 4 * B
HmnR[0, 1, ir] = Delta_0
HmnR[1, 0, ir] = Delta_0
HmnR[2, 3, ir] = Delta_0
HmnR[3, 2, ir] = Delta_0

# 1 0
ir = 1
Irvec[:, ir] = [1, 0, 0]
HmnR[0, 0, ir] = B
HmnR[1, 1, ir] = -B
HmnR[2, 2, ir] = B
HmnR[3, 3, ir] = -B
HmnR[0, 1, ir] =-0.5*zi*A
HmnR[1, 0, ir] =-0.5*zi*A
HmnR[2, 3, ir] = 0.5*zi*A
HmnR[3, 2, ir] = 0.5*zi*A

# 0 1
ir = 2
Irvec[:, ir] = [0, 1, 0]
HmnR[0, 0, ir] = B
HmnR[1, 1, ir] = -B
HmnR[2, 2, ir] = B
HmnR[3, 3, ir] = -B
HmnR[0, 1, ir]=  -A/2
HmnR[1, 0, ir]=   A/2
HmnR[2, 3, ir]=  -A/2
HmnR[3, 2, ir]=   A/2


# -1 0
ir = 3
Irvec[:, ir] = [-1, 0, 0]
HmnR[0, 0, ir] = B
HmnR[1, 1, ir] = -B
HmnR[2, 2, ir] = B
HmnR[3, 3, ir] = -B
HmnR[0, 1, ir]= 0.5*zi*A
HmnR[1, 0, ir]= 0.5*zi*A
HmnR[2, 3, ir]=-0.5*zi*A
HmnR[3, 2, ir]=-0.5*zi*A


# 0 -1
ir = 4
Irvec[:, ir] = [0, -1, 0]
HmnR[0, 0, ir] = B
HmnR[1, 1, ir] = -B
HmnR[2, 2, ir] = B
HmnR[3, 3, ir] = -B
HmnR[0, 1, ir]=   A/2
HmnR[1, 0, ir]=  -A/2
HmnR[2, 3, ir]=   A/2
HmnR[3, 2, ir]=  -A/2

nrpts= ir+1
# Writing to a file
with open('BHZ_hr.dat', 'w') as file:
    file.write('4-band of BHZ model\n')
    file.write('4 !num_wann \n')
    file.write(f'{nrpts} ! nrpts\n')
    file.write(' '.join(f'{x:5d}' for x in ndegen) + '\n')
    for ir in range(nrpts):
        for i in range(4):
            for j in range(4):
                file.write(f"{Irvec[0, ir]:5d}{Irvec[1, ir]:5d}{Irvec[2, ir]:5d}{i+1:5d}{j+1:5d} {HmnR[i, j, ir].real:16.8f} {HmnR[i, j, ir].imag:16.8f}\n")

核验程序读 HR、WCC、左边界和体谱保存列，并读取日志中的自适应积分记录。它检查 Hermiticity、时间反演、反演和有限网格直接隙，再读占据 Kramers 对的四个 TRIM 宇称，按 Fu–Kane 式 (1.1)–(1.2)计算奇偶，与两个 WCC 判定比较。这个宇称核对依赖模型的反演对称性；不要求异质结也具有该对称性。

#!/usr/bin/env python3
"""Read the official BHZ HR and WT outputs; NumPy only; no DFT."""
from pathlib import Path
import itertools, json, re
import numpy as np
root = Path(__file__).resolve().parent
lines = (root / "BHZ_hr.dat").read_text().splitlines()
nw, nr = int(lines[1].split()[0]), int(lines[2].split()[0])
assert (nw, nr) == (4, 5)
# Official generator prints seven unit degeneracies; WT reads five from this line.
deg = np.array([float(x) for x in lines[3].split()][:nr])
raw = np.array([[float(x) for x in s.split()] for s in lines[4:]])
assert raw.shape == (nr*nw*nw, 7) and np.all(deg == 1)
rvec, hr = [], []
for block in raw.reshape(nr, nw*nw, 7):
    assert np.all(block[:, :3] == block[0, :3])
    h = np.zeros((nw,nw), complex)
    for row in block:
        h[int(row[3])-1,int(row[4])-1] = row[5] + 1j*row[6]
    rvec.append(block[0,:3]); hr.append(h)
rvec, hr = np.array(rvec), np.array(hr)
def H(k):
    return np.einsum("r,rij->ij", np.exp(2j*np.pi*(rvec@k))/deg, hr)
T = np.block([[np.zeros((2,2)),np.eye(2)],[-np.eye(2),np.zeros((2,2))]])
P = np.diag([1,-1,1,-1])
gap, herm, tr, inv = float("inf"), 0.0, 0.0, 0.0
for x,y in itertools.product(np.linspace(0,1,101,endpoint=False), repeat=2):
    k = np.array([x,y,0]); h = H(k)
    herm = max(herm, float(np.max(abs(h-h.conj().T))))
    tr = max(tr, float(np.max(abs(T@h.conj()@T.conj().T-H(-k)))))
    inv = max(inv, float(np.max(abs(P@h@P-H(-k)))))
    e = np.linalg.eigvalsh(h); gap = min(gap, float(e[2]-e[1]))
assert herm < 1e-12 and tr < 1e-12 and inv < 1e-12 and gap > 0
parities = []
for k in [[0,0,0],[.5,0,0],[0,.5,0],[.5,.5,0]]:
    e,v = np.linalg.eigh(H(k)); occ = v[:,:2]
    p = np.linalg.eigvalsh(occ.conj().T@P@occ)
    assert np.max(abs(abs(p)-1)) < 1e-12 and abs(p[0]-p[1]) < 1e-12
    parities.append(int(round(p[0])))
z2_parity = int((1-np.prod(parities))//2)
checks = []
for n in [41,81]:
    run = root/f"n{n}-final"; out = (run/"WT.out").read_text()
    assert "ERROR" not in out+(run/"run.out").read_text()
    z2 = int(re.findall(r"Z2 for the plane you choose:\s*(\d+)", out)[-1])
    w = np.loadtxt(run/"wcc.dat")
    assert w.shape == (n,5) and np.isfinite(w).all()
    assert abs(w[0,0]) < 1e-8 and abs(w[-1,0]-.5) < 1e-8
    ends = [float(abs((a[3]-a[4]+.5)%1-.5)) for a in w[[0,-1]]]
    assert max(ends) < 1e-8 and z2 == z2_parity == 1
    l = np.loadtxt(run/"dos.dat_l"); bulk = np.loadtxt(run/"dos.dat_bulk")
    assert l.shape == (n*401,4) and bulk.shape == (n*401,3)
    assert np.isfinite(l).all() and np.isfinite(bulk).all()
    assert np.max(abs(l[:,:2]-bulk[:,:2])) < 1e-10
    # k is an accumulated path length. Middle record corresponds to kx=0.
    center = l.reshape(n,401,4)[n//2]
    j = int(np.argmin(abs(center[:,1])))
    log_l = float(center[j,2]); log_bulk = float(bulk.reshape(n,401,3)[n//2,j,2])
    records = re.findall(r"Wcc integration max_diff, Nk_adaptive\s+([-+0-9.Ee]+)\s+(\d+)", out)
    assert len(records) == n and set(int(x[1]) for x in records) == {64,128}
    counts = {str(k): sum(int(x[1]) == k for x in records) for k in [64,128]}
    tol = float(re.findall(r"wcc_calc_tol\s+([-+0-9.Ee]+)", out)[-1])
    neighbour = float(re.findall(r"wcc_neighbour_tol\s+([-+0-9.Ee]+)", out)[-1])
    max_change = max(float(x[0]) for x in records)
    assert max_change <= tol and tol == .08 and neighbour == .30
    checks.append(dict(mesh=n,input_Nk1=n,input_Nk2=n,
                       wcc_algorithm="dense adaptive integration; initial transverse Nk2 sampling",
                       loop_points_final_counts=counts,loop_max_recorded_change=max_change,
                       wcc_calc_tol=tol,wcc_neighbour_tol=neighbour,
                       z2=z2,wcc_rows=len(w),kramers_endpoint_circular_difference=ends,
                       near_zero_energy_eV=float(center[j,1]),ln_left_ldos=log_l,
                       ln_bulk_ldos=log_bulk,left_over_bulk_ldos=float(np.exp(log_l-log_bulk))))
summary = dict(scope="Official lattice BHZ software check, not a material DFT/Wannier result",
               parameters_model_eV=dict(M=2,B=1,A=1,Delta0=0),occupied_spinor_bands=2,
               full_grid_gap_min_eV=gap,gap_grid="101x101, not a material convergence test",
               hermiticity_max=herm,time_reversal_max=tr,inversion_max=inv,occupied_pair_TRIM_parities=parities,
               z2_inversion_parity=z2_parity,surface_eta_eV=15/401,
               weight_comparison_scope="exp(log-left minus log-bulk) only; projections have different trace dimensions, no equal absolute normalization",
               surface_columns="ln LDOS as written by surfstat.f90; no absolute DOS normalization inferred",
               runs=checks)
(root/"checks.json").write_text(json.dumps(summary,indent=2)+"\n")
print(json.dumps(summary,indent=2))

gnuplot 直接读取实际 WCC 与边界输出；把单段累计路径长度线性映射回输入的 kₓ=−0.5…0.5，保存物理坐标而不依赖原文件的长度单位表题。

本例 KPATH_SLAB 只有这一条直线段，第一列从零累计到 pathmax，所以代码中的 $1/pathmax-0.5 恢复了分数 kₓ；多段路径或曲线不能沿用这个全局缩放。WCC 使用第四、五列的模 1 相位作散点，不补造跨点分支；谱图则把第二列 eV 能量与第三列对数权重直接送入 gnuplot pm3d，保留实际采样点、展宽和固定色标。

# Nk2 labels transverse outputs; dense loop integration is adaptive (64/128 final points).
# WT output columns and periodic path; no recomputation or smoothing.
set encoding utf8
set terminal pngcairo enhanced font "DejaVu Sans,12" size 1400,480
set output 'bhz-wcc-edge.png'
set multiplot layout 1,2 margins 0.08,0.88,0.14,0.91 spacing 0.12
set title 'BHZ model | two occupied spinor bands'
set xlabel 'Transverse k_y (fractional)'
set ylabel 'WCC (mod 1)'
set xrange [0:0.5]
set yrange [0:1]
set xtics 0.1
set ytics 0.25
set key top right
plot 'n81-final/wcc.dat' using 1:4 with points pt 7 ps 0.45 lc rgb '#0072b2' title 'Nk2=81', \
 'n81-final/wcc.dat' using 1:5 with points pt 7 ps 0.45 lc rgb '#0072b2' notitle, \
 'n41-final/wcc.dat' using 1:4 with points pt 6 ps 0.7 lc rgb '#d55e00' title 'Nk2=41', \
 'n41-final/wcc.dat' using 1:5 with points pt 6 ps 0.7 lc rgb '#d55e00' notitle
set title 'Semi-infinite y boundary | eta = 0.0374 eV'
unset key
set xlabel 'Conserved k_x (fractional)'
set ylabel 'Energy (eV)'
set xrange [-0.5:0.5]
set yrange [-2.5:2.5]
set xtics 0.25
set ytics 1
set view map
set palette defined (-5 '#194eff', 0 'white', 5 '#d73027')
set cbrange [-5:5]
set cblabel 'ln LDOS (WT output)'
set pm3d map
stats 'n81-final/dos.dat_l' using 1 nooutput
pathmax = STATS_max
splot 'n81-final/dos.dat_l' using ($1/pathmax-0.5):2:3 with pm3d
unset multiplot