真空参考能级对齐

#!/usr/bin/env python3
"""Export vacuum-referenced band-edge and facing-surface tables from analysis JSON.

Standard library only. Run from this directory after analyze_alignment.py.
"""
from __future__ import annotations

import csv
import json
from pathlib import Path

ROOT = Path(__file__).resolve().parent
SOURCE = ROOT / "alignment-summary.json"


def write_csv(path: Path, columns: list[str], rows: list[dict[str, object]]) -> None:
    with path.open("w", newline="", encoding="utf-8") as handle:
        writer = csv.DictWriter(handle, fieldnames=columns)
        writer.writeheader()
        writer.writerows(rows)


def main() -> None:
    summary = json.loads(SOURCE.read_text(encoding="utf-8"))
    layer_rows: list[dict[str, object]] = []
    for material, layer in summary["layers"].items():
        extrema = {
            "vbm_minus_vacuum_eV": "",
            "cbm_minus_vacuum_eV": "",
        }
        for window in layer["windows"]:
            layer_rows.append({
                "material": material,
                "surface": window["side"],
                "vacuum_window_A": f'{window["lo_A"]:.1f}:{window["hi_A"]:.1f}',
                "vacuum_mean_eV": f'{window["mean_eV"]:.9f}',
                "vacuum_range_eV": f'{window["range_eV"]:.9g}',
                "fermi_minus_vacuum_eV": f'{window["fermi_minus_vacuum_eV"]:.9f}',
                "vbm_minus_vacuum_eV": f'{window.get("vbm_minus_vacuum_eV", ""):.9f}' if "vbm_minus_vacuum_eV" in window else "",
                "cbm_minus_vacuum_eV": f'{window.get("cbm_minus_vacuum_eV", ""):.9f}' if "cbm_minus_vacuum_eV" in window else "",
                "sampled_gap_eV": f'{layer["gap_eV"]:.9f}' if "gap_eV" in layer else "",
                "classification": layer["classification"],
            })
    write_csv(ROOT / "band-edges-vacuum-referenced.csv", [
        "material", "surface", "vacuum_window_A", "vacuum_mean_eV",
        "vacuum_range_eV", "fermi_minus_vacuum_eV",
        "vbm_minus_vacuum_eV", "cbm_minus_vacuum_eV",
        "sampled_gap_eV", "classification",
    ], layer_rows)

    facing = summary["interface_facing_isolated_reference"]
    offsets = [
        {
            "quantity": "CBM(SnSe2, lower-z) - EF(Sr2N, upper-z)",
            "value_eV": f'{facing["cbm_minus_metal_fermi_eV"]:.9f}',
            "meaning": "vacuum-referenced isolated-layer edge offset",
            "scope": facing["scope"],
        },
        {
            "quantity": "EF(Sr2N, upper-z) - VBM(SnSe2, lower-z)",
            "value_eV": f'{facing["metal_fermi_minus_vbm_eV"]:.9f}',
            "meaning": "vacuum-referenced isolated-layer edge offset",
            "scope": facing["scope"],
        },
    ]
    write_csv(ROOT / "facing-surface-offsets.csv", [
        "quantity", "value_eV", "meaning", "scope",
    ], offsets)

    print(f"Wrote {len(layer_rows)} surface rows to band-edges-vacuum-referenced.csv")
    print(f"Wrote {len(offsets)} facing offsets to facing-surface-offsets.csv")


if __name__ == "__main__":
    main()

from __future__ import print_function
import os,re,json,math,hashlib

def sha(p): return hashlib.sha256(open(p,'rb').read()).hexdigest()
def poscar(p):
 l=open(p).read().splitlines(); s=float(l[1]); cell=[[s*float(x) for x in r.split()] for r in l[2:5]]
 if s<=0 or not l[7].lower().startswith('d'): raise ValueError('Expected positive-scale direct-coordinate POSCAR')
 species=l[5].split(); counts=list(map(int,l[6].split())); n=sum(counts)
 xyz=[list(map(float,r.split()[:3])) for r in l[8:8+n]]
 if len(xyz)!=n: raise ValueError('Truncated POSCAR')
 return cell,species,counts,xyz

def geometry():
 ref,sp,n,r=poscar('POSCAR.reference')
 if sp!=['Sn','Se','N','Sr'] or n!=[1,2,1,2]: raise ValueError('Reference species changed')
 out={}
 for name,offset,expect in [('snse2',0,['Sn','Se']),('sr2n',3,['N','Sr'])]:
  cell,species,counts,xyz=poscar(name+'/POSCAR')
  if cell!=ref or species!=expect or counts!=[1,2]: raise ValueError('Cell/species mismatch')
  old=r[offset:offset+3]; shifts=[x[2]-y[2] for x,y in zip(xyz,old)]
  if max(shifts)-min(shifts)>1e-12: raise ValueError('Internal layer geometry changed')
  if any(abs(x[i]-y[i])>1e-12 for x,y in zip(xyz,old) for i in [0,1]): raise ValueError('In-plane registry changed')
  z=[x[2]*cell[2][2] for x in xyz]
  out[name]={'atoms':3,'z_shift_fractional':sum(shifts)/3,'z_min_A':min(z),'z_max_A':max(z),'thickness_A':max(z)-min(z),'empty_height_A':cell[2][2]-(max(z)-min(z))}
  print('%s: 3 atoms; rigid dz=%.12f fractional; thickness=%.9f A; empty height=%.9f A'%(name,out[name]['z_shift_fractional'],out[name]['thickness_A'],out[name]['empty_height_A']))
 if open('snse2/INCAR').read()!=open('sr2n/INCAR').read(): raise ValueError('INCAR protocols differ')
 if open('snse2/KPOINTS').read()!=open('sr2n/KPOINTS').read(): raise ValueError('KPOINTS protocols differ')
 out['common_cell_A']=ref
 print('Common cell, INCAR and KPOINTS: identical')
 return out

def read_eigenval(p):
 f=open(p); head=[f.readline() for _ in range(5)]
 if int(head[0].split()[-1])!=1: raise ValueError('Only ISPIN=1 scalar EIGENVAL is supported')
 ne,nk,nb=map(int,f.readline().split()); points=[]
 for ik in range(nk):
  l=f.readline()
  while l and not l.strip(): l=f.readline()
  k=list(map(float,l.split()))
  if len(k)!=4: raise ValueError('Invalid or truncated k-point header')
  bands=[]
  for ib in range(nb):
   row=list(map(float,f.readline().split()))
   if len(row)!=3 or row[0]!=ib+1: raise ValueError('Invalid band row')
   if not all(not math.isnan(x) and not math.isinf(x) for x in row): raise ValueError('Non-finite eigenvalue')
   if row[2]<-1e-6 or row[2]>1+1e-6: raise ValueError('Unexpected non-spin occupation convention')
   bands.append(row[1:])
  points.append((k,bands))
 f.close()
 w=sum(k[3] for k,rows in points)
 if abs(w-1)>1e-6: raise ValueError('Expected normalized SCF k weights')
 # VASP 5.4.4 scalar EIGENVAL occupations are per spin, in [0,1].
 count=2*sum(k[3]*sum(row[1] for row in rows) for k,rows in points)
 if abs(count-ne)>2e-4: raise ValueError('Weighted occupations do not reproduce NELECT: %.9f vs %d'%(count,ne))
 return ne,nk,nb,points,count

def result(name):
 out=open(name+'/OUTCAR').read()
 if 'aborting loop because EDIFF is reached' not in out or 'General timing and accounting' not in out:
  raise ValueError(name+': SCF convergence/final accounting missing')
 if int(re.search(r'ISPIN\s*=\s*(\d+)',out).group(1))!=1: raise ValueError('Wrong ISPIN')
 ef=float(re.findall(r'E-fermi\s*:\s*([-+0-9.]+)',out)[-1])
 elapsed=float(re.findall(r'Elapsed time \(sec\):\s*([0-9.]+)',out)[-1])
 ne,nk,nb,points,count=read_eigenval(name+'/EIGENVAL')
 if abs(float(re.search(r'NELECT\s*=\s*([0-9.]+)',out).group(1))-ne)>1e-8: raise ValueError('OUTCAR/EIGENVAL electron counts differ')
 crossing=[]; ext=[]; partial=[]
 for ib in range(nb):
  energies=[r[ib][0] for k,r in points]; occ=[r[ib][1] for k,r in points]
  item={'band':ib+1,'emin_eV':min(energies),'emax_eV':max(energies),'occ_min':min(occ),'occ_max':max(occ)}
  ext.append(item)
  if min(energies)<ef<max(energies): crossing.append(ib+1)
  if any(1e-3<x<1-1e-3 for x in occ): partial.append(ib+1)
 res={'fermi_eV':ef,'electrons':ne,'weighted_electrons':count,'nkpoints':nk,'bands':nb,'crossing_bands':crossing,'partially_occupied_bands':partial,'band_extrema':ext,'elapsed_seconds':elapsed,'windows':[],'hashes':{x:sha(name+'/'+x) for x in ['INCAR','POSCAR','KPOINTS','OUTCAR','OSZICAR','EIGENVAL','LOCPOT']}}
 if ne%2==0:
  ib=ne//2
  val=[(r[ib-1][0],i,k[:3]) for i,(k,r) in enumerate(points)]
  con=[(r[ib][0],i,k[:3]) for i,(k,r) in enumerate(points)]
  vbm=max(val);cbm=min(con);gap=cbm[0]-vbm[0]
  if gap>0 and not crossing:
   res.update({'classification':'gapped_on_sampled_mesh','vbm_eV':vbm[0],'cbm_eV':cbm[0],'gap_eV':gap,'vbm_k_fractional':vbm[2],'cbm_k_fractional':cbm[2],'vbm_band':ib,'cbm_band':ib+1})
  else: res['classification']='metallic_on_sampled_mesh' if crossing else 'no_global_gap'
 else:
  if not crossing: raise ValueError(name+': odd-electron scalar case has no sampled crossing; inspect manually')
  res['classification']='metallic_on_sampled_mesh'
 p=json.load(open(name+'/potential-summary.json'))
 if p['source_sha256']!=sha(name+'/LOCPOT'): raise ValueError('Potential summary is stale')
 if len(p['windows'])!=2: raise ValueError('Expected exactly two vacuum windows')
 for side,w in zip(['lower_z','upper_z'],p['windows']):
  if w['range_eV']>0.005: raise ValueError(name+' '+side+': vacuum not flat within 5 meV')
  item=dict(w); item['side']=side; item['fermi_minus_vacuum_eV']=ef-w['mean_eV'];item['workfunction_eV']=w['mean_eV']-ef
  if 'vbm_eV' in res:
   item['vbm_minus_vacuum_eV']=res['vbm_eV']-w['mean_eV'];item['cbm_minus_vacuum_eV']=res['cbm_eV']-w['mean_eV']
  res['windows'].append(item)
 print('%s: NELECT=%d; weighted electrons=%.8f; NKPTS=%d; NBANDS=%d; elapsed=%.3f s'%(name,ne,count,nk,nb,elapsed))
 print('  E_F=%.6f eV; classification=%s; crossing bands=%s'%(ef,res['classification'],crossing))
 for ib in crossing:
  e=ext[ib-1];print('  band %d: Emin=%.6f eV; Emax=%.6f eV; occupation=%.6f:%.6f'%(ib,e['emin_eV'],e['emax_eV'],e['occ_min'],e['occ_max']))
 if 'gap_eV' in res: print('  VBM=%.6f eV; CBM=%.6f eV; sampled gap=%.6f eV'%(res['vbm_eV'],res['cbm_eV'],res['gap_eV']))
 for w in res['windows']:
  print('  %s V_vac=%.9f eV; range=%.6g eV; Phi=%.9f eV'%(w['side'],w['mean_eV'],w['range_eV'],w['workfunction_eV']))
 return res

if __name__=='__main__':
 import sys
 g=geometry()
 if '--geometry-only' in sys.argv:
  json.dump(g,open('geometry-check.json','w'),indent=2,sort_keys=True)
 else:
  results={name:result(name) for name in ['snse2','sr2n']}
  a=results['snse2'];b=results['sr2n']
  if a['classification']!='gapped_on_sampled_mesh' or b['classification']!='metallic_on_sampled_mesh': raise ValueError('The expected semiconductor/metal scope is not supported')
  av=a['windows'][0];bv=b['windows'][1]
  ref={'snse2_side':'lower_z','sr2n_side':'upper_z','cbm_minus_metal_fermi_eV':av['cbm_minus_vacuum_eV']-bv['fermi_minus_vacuum_eV'],'metal_fermi_minus_vbm_eV':bv['fermi_minus_vacuum_eV']-av['vbm_minus_vacuum_eV'],'scope':'Frozen isolated layers, scalar nonmagnetic PBE-D3 (zero damping); not an interface barrier'}
  summary={'geometry':g,'layers':results,'interface_facing_isolated_reference':ref}
  json.dump(summary,open('alignment-summary.json','w'),indent=2,sort_keys=True)
  print('Facing isolated references: CBM(SnSe2)-E_F(Sr2N)=%.6f eV; E_F(Sr2N)-VBM(SnSe2)=%.6f eV'%(ref['cbm_minus_metal_fermi_eV'],ref['metal_fermi_minus_vbm_eV']))

from __future__ import print_function
import json,math
from plane_average import read_grid,cross,dot
summary={}
for name in ['snse2','sr2n']:
 cell,grid,v=read_grid(name+'/CHGCAR')
 volume=abs(dot(cell[0],cross(cell[1],cell[2])))
 nxy=grid[0]*grid[1];height=abs(dot(cell[2],cross(cell[0],cell[1])))/math.sqrt(dot(cross(cell[0],cell[1]),cross(cell[0],cell[1])))
 z=[height*i/grid[2] for i in range(grid[2])]
 density=[sum(v[i*nxy:(i+1)*nxy])/nxy/volume for i in range(grid[2])]
 electrons=sum(v)/len(v)
 expected={'snse2':26,'sr2n':25}[name]
 if abs(electrons-expected)>1e-4: raise ValueError('CHGCAR integral does not match NELECT')
 row={'electrons':electrons,'volume_A3':volume,'grid':grid,'windows':[]}
 print('%s: integrated valence electrons=%.9f; volume=%.9f A^3'%(name,electrons,volume))
 for lo,hi in [(6.,10.),(29.,33.)]:
  a=[n for zz,n in zip(z,density) if lo<=zz<=hi]
  r={'lo_A':lo,'hi_A':hi,'mean_e_per_A3':sum(a)/len(a),'max_abs_e_per_A3':max(abs(x) for x in a)}
  row['windows'].append(r)
  print('  z=%.1f:%.1f A: mean density=%.6g; max abs density=%.6g e/A^3'%(lo,hi,r['mean_e_per_A3'],r['max_abs_e_per_A3']))
 summary[name]=row
json.dump(summary,open('vacuum-density.json','w'),indent=2,sort_keys=True)

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)

[bcgong@localhost band_alignment]$ python split_paw.py
[bcgong@localhost band_alignment]$ cat snse2/POTCAR.identity.txt sr2n/POTCAR.identity.txt
TITEL  = PAW_PBE Sn_d 06Sep2000
POMASS =  118.710; ZVAL   =   14.000    mass and valenz
TITEL  = PAW_PBE Se 06Sep2000
POMASS =   78.960; ZVAL   =    6.000    mass and valenz
TITEL  = PAW_PBE N 08Apr2002
POMASS =   14.001; ZVAL   =    5.000    mass and valenz
TITEL  = PAW_PBE Sr_sv 07Sep2000
POMASS =   87.620; ZVAL   =   10.000    mass and valenz

watch -n 5 'squeue -o "%.18i %.14j %.10u %.2t %.6C %.10M"'
tail -f out

[bcgong@localhost band_alignment]$ grep -A 5 "FREE ENERGIE" snse2/OUTCAR
FREE ENERGIE OF THE ION-ELECTRON SYSTEM (eV)
  ---------------------------------------------------
  free  energy   TOTEN  =       -12.02561657 eV

  energy  without entropy=      -12.02561626  energy(sigma->0) =      -12.02561641

[bcgong@localhost band_alignment]$ ls -lh snse2/{INCAR,POSCAR,KPOINTS,OUTCAR,OSZICAR,EIGENVAL,LOCPOT,CHGCAR,WAVECAR}
-rw-rw-r-- 1 bcgong bcgong  37M Sep 22 23:33 snse2/CHGCAR
-rw-rw-r-- 1 bcgong bcgong  35K Sep 22 23:33 snse2/EIGENVAL
-rw-rw-r-- 1 bcgong bcgong  369 Sep 22 23:30 snse2/INCAR
-rw-rw-r-- 1 bcgong bcgong   57 Sep 22 23:30 snse2/KPOINTS
-rw-rw-r-- 1 bcgong bcgong  37M Sep 22 23:33 snse2/LOCPOT
-rw-rw-r-- 1 bcgong bcgong 2.4K Sep 22 23:33 snse2/OSZICAR
-rw-rw-r-- 1 bcgong bcgong 135K Sep 22 23:33 snse2/OUTCAR
-rw-rw-r-- 1 bcgong bcgong  464 Sep 22 23:30 snse2/POSCAR
-rw-rw-r-- 1 bcgong bcgong    0 Sep 22 23:31 snse2/WAVECAR

明确的共同晶胞与六原子结构
  └─ 分出两层 → 保留内部形变与表面方向 → 各自居中
       └─ 相同协议的两份独立静态 SCF
            ├─ OUTCAR / OSZICAR：电子收敛与结束验收
            ├─ EIGENVAL：半导体带边或金属费米面交叉
            └─ LOCPOT + CHGCAR：双侧平坦真空区间
                 └─ 每个能级减去对应表面的真空势
                      └─ 孤立层参考表 → 后续直接检查界面体系

# Actual frozen-layer SCF table; no fitting and no synthetic metal band edges.
# Download band-edges-vacuum-referenced.csv beside this script, then run gnuplot.
if (!exists('datafile')) datafile = 'band-edges-vacuum-referenced.csv'
if (!exists('prefix')) prefix = 'frozen-facing-alignment'
set datafile separator comma
stats datafile using (strcol(1) eq 'snse2' && strcol(2) eq 'lower_z' ? column(7) : 1/0) nooutput
if (STATS_records != 1) { print 'Need exactly one SnSe2 lower-z VBM'; exit error }
vbm = STATS_mean
stats datafile using (strcol(1) eq 'snse2' && strcol(2) eq 'lower_z' ? column(8) : 1/0) nooutput
if (STATS_records != 1) { print 'Need exactly one SnSe2 lower-z CBM'; exit error }
cbm = STATS_mean
stats datafile using (strcol(1) eq 'sr2n' && strcol(2) eq 'upper_z' ? column(6) : 1/0) nooutput
if (STATS_records != 1) { print 'Need exactly one Sr2N upper-z EF'; exit error }
ef = STATS_mean
offset = cbm - ef
set encoding utf8
set border 2 lc rgb '#444444'
set tics nomirror out
set xrange [0.5:2.8]
set yrange [-6.65:0.65]
set ylabel 'Energy relative to the selected vacuum (eV)'
set xtics ('SnSe₂  lower-z' 1, 'Sr₂N  upper-z' 2.2) scale 0
set ytics 1
set key off
set grid ytics lc rgb '#dddddd'
set object 1 rect from 0.75,vbm to 1.25,cbm fc rgb '#e9f0f5' fs solid 1 noborder behind
set arrow 1 from 0.75,vbm to 1.25,vbm nohead lw 3 lc rgb '#205a83'
set arrow 2 from 0.75,cbm to 1.25,cbm nohead lw 3 lc rgb '#205a83'
set arrow 3 from 1.95,ef to 2.45,ef nohead lw 3 lc rgb '#ae542f'
set arrow 4 from 1.60,ef to 1.60,cbm heads size screen 0.012,15 lw 1.3 lc rgb '#555555'
set arrow 5 from 1.27,cbm to 1.60,cbm nohead dt 2 lc rgb '#777777'
set arrow 6 from 1.60,ef to 1.93,ef nohead dt 2 lc rgb '#777777'
set label 1 sprintf('VBM  %.6f',vbm) at 0.78,vbm-0.24 left font ',12'
set label 2 sprintf('CBM  %.6f',cbm) at 0.78,cbm+0.25 left font ',12'
set label 3 sprintf('EF  %.6f',ef) at 2.2,ef+0.25 center font ',12'
set label 4 'CBM − EF' at 1.69,(ef+cbm)/2+0.12 left font ',12'
set label 6 sprintf('= %.6f eV',offset) at 1.69,(ef+cbm)/2-0.18 left font ',12'
set label 5 'Vacuum = 0' at 0.58,0.18 left font ',12'
set title 'Frozen isolated layers before contact' font ',17'
do for [export_index=1:3] {
 if (export_index==1) { set terminal svg size 960,680 font 'Liberation Sans,16'; set output prefix.'.svg' }
 if (export_index==2) { set terminal pdfcairo enhanced color size 7in,5in font 'Liberation Sans,11'; set output prefix.'.pdf' }
 if (export_index==3) { set terminal pngcairo size 960,680 font 'Liberation Sans,16'; set output prefix.'.png' }
 plot 0 with lines lw 1.3 dt 2 lc rgb '#777777' notitle
 unset output
}
print sprintf('VBM=%.9f; CBM=%.9f; Sr2N EF=%.9f; CBM-EF=%.9f eV',vbm,cbm,ef,offset)