态密度 DOS / PDOS

from atlas_plot_style import install as install_atlas_style
install_atlas_style()
from pathlib import Path
import argparse,csv
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
ROOT=Path(__file__).resolve().parent
OUT=ROOT/'plots';OUT.mkdir(exist_ok=True)
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'figure.dpi':130,'savefig.dpi':300,'axes.titleweight':'bold'})
BLUE='#0072b2';ORANGE='#e69f00';GRAY='#667080'

def read(name):return np.genfromtxt(ROOT/name,delimiter=',',names=True,encoding='utf-8')
def numeric_rows(name):
    with (ROOT/name).open(newline='') as f:
        return [{k:(v if k=='directory' else float(v)) for k,v in row.items()}
                for row in csv.DictReader(f)]
def save(fig,name):
    fig.savefig(OUT/(name+'.png'),bbox_inches='tight');fig.savefig(OUT/(name+'.svg'),bbox_inches='tight');plt.close(fig);print('plots/'+name+'.png')
def clean(ax):ax.grid(alpha=.18);ax.set_axisbelow(True)
def convergence():
    rows=list(csv.DictReader((ROOT/'convergence.csv').open()))
    fig,axes=plt.subplots(1,3,figsize=(12,3.6),layout='constrained')
    for ax,param,label in zip(axes,['ecutwfc','ecutrho','kmesh'],['Wavefunction cutoff (Ry)','Charge-density cutoff (Ry)','Uniform mesh n × n × n']):
        r=[x for x in rows if x['parameter']==param];x=[float(v['setting']) for v in r];y=[float(v['delta_meV_per_atom_vs_last']) for v in r]
        ax.plot(x,y,'o-',color=BLUE);ax.set_xlabel(label);ax.set_ylabel('ΔE vs last point (meV/atom)');clean(ax)
        if param=='kmesh':ax.set_yscale('symlog',linthresh=.1);ax.axhline(1,color=ORANGE,ls='--',label='Example comparison: 1 meV/atom');ax.legend(fontsize=8)
        ax.set_title({'ecutwfc':'Fixed 640 Ry, 8³ mesh','ecutrho':'Fixed 60 Ry, 8³ mesh','kmesh':'Fixed 60 / 640 Ry'}[param],fontsize=11)
    save(fig,'convergence')
def relax():
    r=read('relax/relaxation.csv')
    fig,axes=plt.subplots(1,2,figsize=(9.6,3.6),layout='constrained')
    e_diff = (r['energy_Ry'] - r['energy_Ry'][-1]) * 13605.693122994
    axes[0].plot(r['scf_cycle'], e_diff, 'o-', color=BLUE, lw=1.5, ms=5, zorder=3)
    axes[0].fill_between(r['scf_cycle'], 0, e_diff, color=BLUE, alpha=0.12, zorder=1)
    axes[0].set_ylabel('Energy above final geometry (meV/cell)')
    axes[0].set_title('Si: Monotonic Energy Descent', fontweight='bold', fontsize=10.5)
    for cycle, de in zip(r['scf_cycle'][:-1], e_diff[:-1]):
        axes[0].annotate(f"{de:.1f} meV", xy=(cycle, de), xytext=(0, 6),
                         textcoords='offset points', ha='center', fontsize=7.5, color=BLUE)
    axes[0].annotate("Final E₀ = -22.838592 Ry", xy=(r['scf_cycle'][-1], 0), xytext=(-55, 18),
                     textcoords='offset points', fontsize=7.5, color='#162232',
                     arrowprops=dict(arrowstyle='->', color='#718096', lw=0.6))
    
    # Total force plot with log-scale threshold lines
    f_vals = np.array(r['total_force_Ry_per_Bohr'])
    # For display of the 0.0 final step, set a small visual floor for the arrow
    f_disp = np.maximum(f_vals, 1e-5)
    axes[1].plot(r['scf_cycle'][:-1], f_vals[:-1], 's-', color=ORANGE, lw=1.5, ms=5, label='Total Force', zorder=3)
    axes[1].plot([r['scf_cycle'][-2], r['scf_cycle'][-1]], [f_vals[-2], f_disp[-1]], 's--', color=ORANGE, lw=1.0, zorder=2)
    axes[1].scatter([r['scf_cycle'][-1]], [f_disp[-1]], marker='o', facecolors='white', edgecolors=ORANGE, s=40, zorder=4)
    axes[1].axhline(1e-3, color='#c0392b', ls='--', lw=0.9, label=r'QE Default $10^{-3}$ Ry/Bohr', zorder=1)
    axes[1].axhline(1e-4, color='#27ae60', ls=':', lw=0.9, label=r'Stringent $10^{-4}$ Ry/Bohr', zorder=1)
    axes[1].set_yscale('log')
    axes[1].set_ylim(5e-6, 1.2e-1)
    axes[1].set_ylabel('Total Force (Ry/Bohr, log scale)')
    axes[1].set_title('Si: BFGS Force Convergence', fontweight='bold', fontsize=10.5)
    axes[1].annotate('Printed as 0.000000\n(Force strictly zero by symmetry)',
                     xy=(r['scf_cycle'][-1], f_disp[-1]), xytext=(-120, 28),
                     textcoords='offset points', fontsize=7.5, color='#8c3b00',
                     arrowprops=dict(arrowstyle='->', color=ORANGE, lw=0.7))
    axes[1].legend(frameon=True, facecolor='#f8fafc', edgecolor='#e2e8f0', fontsize=7.5, loc='upper right')
    
    for label, ax in zip('ab', axes):
        ax.set_xlabel('BFGS Step / SCF Cycle')
        ax.set_xticks(r['scf_cycle'])
        clean(ax)
        ax.text(-0.14, 1.05, label, transform=ax.transAxes, fontweight='bold', fontsize=10)
    save(fig, 'relax')

def gap():
    r=numeric_rows('gap-results.csv');base=[x for x in r if x['directory'] in ['gap12','gap18-cg','gap24-cg']]
    fig,axes=plt.subplots(1,2,figsize=(10,3.6),layout='constrained')
    axes[0].plot([12,18,24],[x['gap_eV'] for x in base],'o-',color=BLUE);axes[0].set_xticks([12,18,24]);axes[0].set_xlabel('Uniform NSCF mesh n × n × n');axes[0].set_ylabel('Sampled indirect gap (eV)');axes[0].set_title('The sampled minimum need not vary monotonically',fontsize=10)
    for n,x in zip([12,18,24],base):axes[0].annotate(f"{max(x[f'cbm_k{axis}_tpiba'] for axis in 'xyz'):.4f} × 2π/a",(n,x['gap_eV']),xytext=(0,7),textcoords='offset points',ha='center',fontsize=8)
    m=read('mass/longitudinal.csv');fit=numeric_rows('mass/mass-fits.csv')[1]
    dense=next(x for x in r if x['directory']=='gap24-k12-cg');axes[1].plot(m['kx_tpiba'],m['band5_eV']-dense['vbm_eV'],color=BLUE)
    axes[1].scatter([fit['minimum_k_tpiba']],[fit['minimum_energy_eV']-dense['vbm_eV']],color=ORANGE,zorder=4)
    axes[0].margins(y=.20)
    axes[1].set_xlabel('kx (2π/a), ky = kz = 0');axes[1].set_ylabel('Conduction energy − VBM (eV)');axes[1].set_title('Local Γ–X valley refinement; 12³ parent density',fontsize=10)
    for ax in axes:clean(ax)
    save(fig,'band-gap')
def mass():
    r=read('mass/longitudinal.csv');fits=numeric_rows('mass/mass-fits.csv');fit=fits[1]
    fig,axes=plt.subplots(1,2,figsize=(10,3.6),layout='constrained')
    dk=r['kx_inv_A']-fit['minimum_k_inv_A'];take=np.abs(dk)<.04
    axes[0].scatter(dk[take],(r['band5_eV'][take]-fit['minimum_energy_eV'])*1000,s=24,color=BLUE,label='Actual QE eigenvalues')
    q=np.linspace(-.035,.035,101);axes[0].plot(q,fit['quadratic_A_eVA2']*q*q*1000,color=ORANGE,label='Quadratic fit, ±0.02 Å⁻¹')
    axes[0].set_xlabel('kx − valley minimum (Å⁻¹)');axes[0].set_ylabel('Energy above minimum (meV)');axes[0].legend(fontsize=8)
    axes[1].plot([f['half_window_inv_A'] for f in fits],[f['mass_over_me'] for f in fits],'o-',color=BLUE)
    axes[1].set_xlabel('Fit half-window (Å⁻¹)');axes[1].set_ylabel('Longitudinal electron mass / me');axes[1].set_xticks([.01,.02,.03]);axes[1].ticklabel_format(useOffset=False,axis='y')
    for ax in axes:clean(ax)
    save(fig,'effective-mass')
def band3d():
    r=read('band3d/cube.csv');emin=r['band5_eV'].min();fig=plt.figure(figsize=(11,4.4),layout='constrained')
    for i,(key,value,xkey,ykey,title) in enumerate([('kz_tpiba',0.,'kx_tpiba','ky_tpiba','kz = 0 plane'),('kx_tpiba',.85,'ky_tpiba','kz_tpiba','kx = 0.85 × 2π/a plane')]):
        d=r[np.isclose(r[key],value)];xs=np.unique(d[xkey]);ys=np.unique(d[ykey]);X,Y=np.meshgrid(xs,ys,indexing='ij');Z=np.empty_like(X)
        for u,x in enumerate(xs):
            for v,y in enumerate(ys):Z[u,v]=(d['band5_eV'][np.isclose(d[xkey],x)&np.isclose(d[ykey],y)][0]-emin)*1000
        ax=fig.add_subplot(1,2,i+1,projection='3d');ax.plot_surface(X,Y,Z,cmap='viridis',linewidth=.2,edgecolor=(0,0,0,.15),antialiased=True)
        ax.set_xlabel(xkey.replace('_tpiba','')+' (2π/a)');ax.set_ylabel(ykey.replace('_tpiba','')+' (2π/a)');ax.set_zlabel('');ax.text2D(0.02, 0.91, 'E − sampled minimum (meV)', transform=ax.transAxes, fontsize=9);ax.set_title(title);ax.view_init(elev=25,azim=-125);ax.xaxis.set_major_locator(plt.MaxNLocator(4));ax.yaxis.set_major_locator(plt.MaxNLocator(4));ax.zaxis.set_major_locator(plt.MaxNLocator(4))
    fig.suptitle('Two cuts through an actual 11 × 9 × 9 local k-point cube',fontsize=12)
    save(fig,'band-3d')
def population():
    r=read('population-cg/lowdin.csv');fig,ax=plt.subplots(figsize=(6,4),layout='constrained');x=np.arange(len(r))
    ax.bar(x,r['s_electrons'],color=BLUE,label='s');ax.bar(x,r['p_electrons'],bottom=r['s_electrons'],color=ORANGE,label='p')
    for i,q in enumerate(r['total_electrons']):ax.text(i,q+.03,f'{q:.4f}',ha='center')
    ax.axhline(4,color=GRAY,ls='--',lw=1,label='4 valence electrons / Si');ax.set_xticks(x,['Si 1','Si 2']);ax.set_ylabel('Projected electrons');ax.set_ylim(0,4.5);ax.legend(ncol=3,fontsize=9,loc='lower center',bbox_to_anchor=(.5,1.));ax.set_title('Equivalent atoms; spilling = 0.0092',pad=38);save(fig,'population-analysis')
def fatband():
    r=read('bands-cg/fatband.csv');vbm=r['energy_eV'][r['iband']==4].max();fig,axes=plt.subplots(1,2,figsize=(12,4.6),sharey=True,layout='constrained')
    idx=[1,25,37,49,73,97,121];ticks=[r['path_distance_tpiba'][r['ik']==i][0] for i in idx];labels=['Γ','X','W','K','Γ','L','X']
    for ax,weight,color,title in zip(axes,['Si_s_weight','Si_p_weight'],[BLUE,ORANGE],['Si s projection','Si p projection']):
        for ib in range(1,9):
            q=r[r['iband']==ib];ax.plot(q['path_distance_tpiba'],q['energy_eV']-vbm,color='#b6bdc7',lw=.65,zorder=1);ax.scatter(q['path_distance_tpiba'],q['energy_eV']-vbm,s=34*q[weight],color=color,alpha=.75,edgecolors='none',zorder=2)
        for tick in ticks:ax.axvline(tick,color='#e0e3e8',lw=.65,zorder=0)
        ax.axhline(0,color=GRAY,ls='--',lw=.7);ax.set_xticks(ticks,labels);ax.set_xlim(ticks[0],ticks[-1]);ax.set_ylim(-13,7);ax.set_title(title);ax.set_ylabel('Energy − VBM (eV)')
        handles=[ax.scatter([],[],s=34*w,color='#555555',edgecolors='none',label=f'{w:g}') for w in (.25,.5,1.)]
        ax.legend(handles=handles,title='Projection weight (point area)',ncol=3,
                  loc='lower left',bbox_to_anchor=(0,1.015),frameon=False,fontsize=8)
        ax.set_title('',loc='center')
        ax.set_title(title,loc='left',pad=55)
    save(fig,'fatband')
def bands():
    r=read('bands-cg/fatband.csv');vbm=r['energy_eV'][r['iband']==4].max();fig,ax=plt.subplots(figsize=(8,4.6),layout='constrained');idx=[1,25,37,49,73,97,121];ticks=[r['path_distance_tpiba'][r['ik']==i][0] for i in idx]
    for ib in range(1,9):
        q=r[r['iband']==ib];ax.plot(q['path_distance_tpiba'],q['energy_eV']-vbm,color=BLUE,lw=1)
    for tick in ticks:ax.axvline(tick,color='#dce1e8',lw=.7)
    ax.axhline(0,color=GRAY,ls='--',lw=.7);ax.set_xticks(ticks,['Γ','X','W','K','Γ','L','X']);ax.set_xlim(ticks[0],ticks[-1]);ax.set_ylim(-13,7);ax.set_ylabel('Energy − VBM (eV)');ax.set_title('Si, PBE, fixed example cell; no SOC');save(fig,'bands')
def dos():
    r=np.loadtxt(ROOT/'dos-cg/si.dos.dat');g=numeric_rows('gap-results.csv');vbm=next(x['vbm_eV'] for x in g if x['directory']=='gap24-cg')
    fig,ax=plt.subplots(figsize=(7.5,3.7),layout='constrained');ax.plot(r[:,0]-vbm,r[:,1],color=BLUE);ax.fill_between(r[:,0]-vbm,r[:,1],alpha=.15,color=BLUE);ax.axvline(0,color=GRAY,ls='--');ax.set_xlim(-13,8);ax.set_xlabel('Energy − VBM (eV)');ax.set_ylabel('DOS (states/eV/cell)');ax.set_title('24³ NSCF mesh; Gaussian broadening 0.01 Ry');clean(ax);save(fig,'dos')

if __name__=='__main__':
    parser=argparse.ArgumentParser();parser.add_argument('figure',choices=['convergence','relax','gap','mass','band3d','population','fatband','bands','dos','all']);args=parser.parse_args()
    funcs={'convergence':convergence,'relax':relax,'gap':gap,'mass':mass,'band3d':band3d,'population':population,'fatband':fatband,'bands':bands,'dos':dos}
    for name,fun in funcs.items():
        if args.figure in [name,'all']:fun()

&CONTROL
  calculation = 'scf'
  restart_mode = 'from_scratch'
  prefix = 'sc2c'
  outdir = './out/'
  pseudo_dir = '../pseudos'
  verbosity = 'high'
  tprnfor = .true.
  tstress = .true.
/
&SYSTEM
  ibrav = 0
  nat = 3
  ntyp = 2
  ecutwfc = 100
  ecutrho = 800
  input_dft = 'vdw-DF3-opt1'
  occupations = 'smearing'
  smearing = 'gaussian'
  degauss = 3.7d-3
  tot_charge = 0.0
  nr1 = 60
  nr2 = 60
  nr3 = 720
  nr1s = 45
  nr2s = 45
  nr3s = 486
/
&ELECTRONS
  conv_thr = 1.0d-12
  mixing_beta = 0.4
/
ATOMIC_SPECIES
Sc  44.956  Sc.pbe-spn-kjpaw_psl.1.0.0.UPF
C  12.011  C.pbe-n-kjpaw_psl.1.0.0.UPF
CELL_PARAMETERS (angstrom)
3.358477221000 -0.000000000000 0.000000000000
-1.679238611000 2.908526592000 0.000000000000
0.000000000000 0.000000000000 40.000000000000
ATOMIC_POSITIONS (angstrom)
C 0.000000000000 0.000000000000 17.688482212000
Sc -0.000000000501 1.939017728097 16.567204708000
Sc 1.679238610501 0.969508863903 18.892666444000
K_POINTS automatic
16 16 1 0 0 0

"""Export one scalar, nonmagnetic QE DOS case from its actual XML and files."""
from pathlib import Path
import argparse
import csv
import math
import re
import xml.etree.ElementTree as ET

HA_EV = 27.211386245988
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--xml', type=Path, required=True)
parser.add_argument('--dos-input', type=Path, required=True)
parser.add_argument('--dos-data', type=Path, required=True)
parser.add_argument('--output-dir', type=Path, required=True)
args = parser.parse_args()
root = ET.parse(args.xml).getroot()
for element in root.iter():
    element.tag = element.tag.split('}')[-1]
band_structure = root.find('./output/band_structure')
if band_structure is None:
    raise ValueError('Missing QE output band_structure')
for flag in ['lsda', 'noncolin', 'spinorbit']:
    element = band_structure.find(flag)
    if element is not None and element.text.strip().lower() in ['true', 't', '.true.']:
        raise ValueError('This scalar NM example does not implement spinor/LSDA export')
ef = float(band_structure.find('fermi_energy').text) * HA_EV
nelec = float(band_structure.find('nelec').text)
nbnd = int(band_structure.find('nbnd').text)
points = band_structure.findall('ks_energies')
if len(points) != int(band_structure.find('nks').text):
    raise ValueError('XML k-point count mismatch')

eigen_rows = []
weights = []
weighted_occupations = 0.0
for index, point in enumerate(points, start=1):
    kpoint = point.find('k_point')
    coordinates = [float(value) for value in kpoint.text.split()]
    weight = float(kpoint.attrib['weight'])
    energies = [float(value) * HA_EV for value in point.find('eigenvalues').text.split()]
    occupations = [float(value) for value in point.find('occupations').text.split()]
    if len(coordinates) != 3 or len(energies) != nbnd or len(occupations) != nbnd:
        raise ValueError('XML band/coordinate count mismatch')
    if weight <= 0 or not all(math.isfinite(v) for v in coordinates + energies + occupations + [weight]):
        raise ValueError('Invalid k-point data')
    weights.append(weight)
    weighted_occupations += weight * sum(occupations)
    for band, energy in enumerate(energies, start=1):
        eigen_rows.append([index, *coordinates, weight, band, energy, energy - ef])
if abs(sum(weights) - 2.0) > 1e-8 or abs(weighted_occupations - nelec) > 1e-6:
    raise ValueError('Scalar NM weight/electron conservation check failed')

namelist = args.dos_input.read_text()
def number(name):
    match = re.search(r'\b' + name + r'\s*=\s*([+\-\d.eEdD]+)', namelist, re.I)
    if match is None:
        raise ValueError('Missing DOS input field: ' + name)
    return float(match.group(1).replace('D', 'e').replace('d', 'e'))
emin, emax, step = number('Emin'), number('Emax'), number('DeltaE')
if step <= 0 or emax <= emin:
    raise ValueError('Invalid DOS energy window')
raw = []
for line in args.dos_data.read_text().splitlines():
    if not line.strip() or line.lstrip().startswith('#'):
        continue
    values = [float(value.replace('D', 'e')) for value in line.split()]
    if len(values) != 3 or not all(math.isfinite(v) for v in values):
        raise ValueError('Invalid native DOS row')
    raw.append(values)
if len(raw) != int(round((emax - emin) / step)) + 1:
    raise ValueError('Native DOS row count does not match this example window')
dos_rows = []
for index, (printed_energy, dos, cumulative) in enumerate(raw):
    exact_energy = emin + index * step
    if abs(printed_energy - exact_energy) > 0.00050001:
        raise ValueError('Native F8.3 energy serialization mismatch')
    dos_rows.append([printed_energy, printed_energy - ef, exact_energy,
                     exact_energy - ef, dos, cumulative])

args.output_dir.mkdir(parents=True, exist_ok=True)
def write_csv(name, fields, rows):
    with (args.output_dir / name).open('w', newline='') as stream:
        writer = csv.writer(stream)
        writer.writerow(fields)
        writer.writerows(rows)
write_csv('dos_shifted.csv', ['native_printed_energy_eV_F8p3',
          'native_printed_energy_minus_own_EF_eV', 'exact_input_grid_energy_eV',
          'exact_input_grid_energy_minus_own_EF_eV',
          'DOS_states_per_eV_per_cell_spin_summed', 'native_window_cumulative_DOS_states'], dos_rows)
write_csv('eigenvalues_k16.csv', ['irreducible_k_index', 'kx_2pi_over_alat',
          'ky_2pi_over_alat', 'kz_2pi_over_alat', 'raw_QE_weight_including_spin',
          'band_index', 'energy_eV_raw', 'energy_minus_own_EF_eV'], eigen_rows)
print(f"Fermi energy: {ef:.12f} eV")
print(f"Electrons: {nelec:g}; bands: {nbnd}; irreducible k points: {len(points)}")
print(f"Spin-inclusive weight sum: {sum(weights):.12f}; weighted occupations: {weighted_occupations:.12f}")
print(f"Wrote {len(dos_rows)} DOS rows and {len(eigen_rows)} eigenvalue rows to {args.output_dir}")
print("DOS and window cumulative values retained; energies shifted by the XML Fermi energy.")

# Run from sc2c-dos/: gnuplot source/plot_dos.gp
# Col4: exact input energy minus XML EF, eV. Col5: spin-summed DOS, states/eV/cell.
# All 101 rows are plotted as straight segments; no smoothing, rescaling or interpolation.
set datafile separator comma
set encoding utf8
set border 3
set tics nomirror out
set xrange [0:*]
set yrange [-0.5:0.5]
set xlabel "DOS (states / eV / cell, both spins)"
set ylabel "E - E_F (eV)"
set title "Sc_2C: fixed three-atom geometry, NM / no SOC"
set key off
set arrow 1 from graph 0, first 0 to graph 1, first 0 nohead dt 2 lc rgb "#666666"
set label 1 "D(E_F) = 5.433" at first 5.433, first 0 offset 1,1 font ",10"
set terminal svg size 660,700 enhanced font "DejaVu Sans,13"
set output "figures/sc2c-dos.svg"
plot "data/dos_shifted.csv" every ::1 using 5:4 with lines lw 2 lc rgb "#215d96", \
     "data/dos_shifted.csv" every ::1 using (abs($4)<1e-8?$5:1/0):4 with points pt 7 ps 1 lc rgb "#a6422b"
set terminal pdfcairo enhanced color size 5.5in,5.8in font "DejaVu Sans,11"
set output "figures/sc2c-dos.pdf"
replot
unset output

处理逻辑是保留前表的三种单位和四个状态名称，只对Sc₂C这一相同层定义作配对。输入是已经由原长表提取的 frozen-window.csv，不把图像上的高度再读成数字；脚本用全精度CSV值计算百分比，并拒绝覆盖源或已有结果。可把下面的具体需求交给编程助手：

读取既有 frozen-window.csv，要求四个唯一冻结状态及有限、非负的Sc2C读数。
在0%和+1.5%下分别比较孤立层→界面内层，再比较界面内层0%→+1.5%。
对D(EF)与[-0.1,0]窗口谱重分别保留参照、目标、差值和100*(target/reference-1)。
D的原单位为states/eV/cell，窗口谱重为states/cell，百分比不是电子转移量。
仅用Python标准库输出新CSV和终端记录，不改源CSV，不画重复的柱图，不调用DFT。

完整比较源码 · 全精度结果 · 实际终端记录。在脚本和原CSV所在工作目录运行：

python3 compare_frozen_spectra.py frozen-window.csv --output frozen-comparisons.csv

实际输出：

Contact at 0%: D(EF) change=-39.1251% ; window change=-58.6958%
Contact at +1.5%: D(EF) change=-54.0299% ; window change=-65.8264%
Heterostructure strain 0% to +1.5%: D(EF) change=-21.1099% ; window change=+5.0139%
These are changes of broadened projected spectra, not transferred electrons.

#!/usr/bin/env python3
"""Compare existing frozen Sc2C spectral readouts; no DFT or source writes."""
from pathlib import Path
import argparse
import csv
import math

D = "sc2c_pdos_at_EF_states_per_eV_cell"
W = "sc2c_pdos_window_weight_states_per_cell"
PAIRS = [
    ("Contact at 0%", "Isolated Sc2C 0%", "Heterostructure 0%"),
    ("Contact at +1.5%", "Isolated Sc2C +1.5%", "Heterostructure +1.5%"),
    ("Heterostructure strain 0% to +1.5%", "Heterostructure 0%", "Heterostructure +1.5%"),
]

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("source", type=Path)
    parser.add_argument("--output", type=Path, required=True)
    args = parser.parse_args()
    if args.output.resolve() == args.source.resolve():
        raise ValueError("Output must differ from the source")
    with args.source.open(newline="") as handle:
        items = list(csv.DictReader(handle))
    rows = {item["state"]: item for item in items}
    if len(items) != 4 or len(rows) != 4:
        raise ValueError("Expected four unique frozen states")
    results = []
    for label, ref, target in PAIRS:
        record = dict(comparison=label, reference=ref, target=target)
        for key, name in [(D, "D_EF"), (W, "window_weight")]:
            x, y = float(rows[ref][key]), float(rows[target][key])
            if not all(math.isfinite(v) for v in [x, y]) or x <= 0 or y < 0:
                raise ValueError("Invalid recorded spectral quantity")
            record[name + "_reference"] = x
            record[name + "_target"] = y
            record[name + "_difference"] = y - x
            record[name + "_change_percent"] = 100 * (y / x - 1)
        results.append(record)
    args.output.parent.mkdir(parents=True, exist_ok=True)
    with args.output.open("x", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=list(results[0]))
        writer.writeheader()
        writer.writerows(results)
    for row in results:
        print("{}: D(EF) change={:+.4f}% ; window change={:+.4f}%".format(
            row["comparison"], row["D_EF_change_percent"],
            row["window_weight_change_percent"]))
    print("These are changes of broadened projected spectra, not transferred electrons.")

if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Read existing ZrCl2/Sc2C tables. No DFT jobs or source writes."""
from pathlib import Path
import argparse, csv, json, math, re


def interp(rows, key, x):
    for a, b in zip(rows, rows[1:]):
        xa, xb = float(a['energy_minus_fermi_eV']), float(b['energy_minus_fermi_eV'])
        if xa <= x <= xb:
            f = (x-xa)/(xb-xa)
            return float(a[key])*(1-f)+float(b[key])*f
    raise ValueError('Energy outside recorded grid')


def integrate(rows, key, lo=-0.1, hi=0.0):
    pts = [(lo, interp(rows,key,lo))]
    pts += [(float(r['energy_minus_fermi_eV']),float(r[key])) for r in rows
            if lo < float(r['energy_minus_fermi_eV']) < hi]
    pts += [(hi, interp(rows,key,hi))]
    return sum((b[0]-a[0])*(a[1]+b[1])/2 for a,b in zip(pts,pts[1:]))


def read_frozen(root):
    groups = {}
    with (root/'frozen-controls-pdos_20260929/frozen_pdos_long.csv').open() as f:
        for r in csv.DictReader(f):
            if not all(math.isfinite(float(r[k])) for k in
                       ['energy_minus_fermi_eV','sc2c_pdos','zrcl2_pdos','total_projected_pdos']):
                raise ValueError('Nonfinite PDOS')
            groups.setdefault(r['state'],[]).append(r)
    out=[]
    for state, rows in groups.items():
        xs=[float(r['energy_minus_fermi_eV']) for r in rows]
        if any(b <= a for a,b in zip(xs,xs[1:])): raise ValueError('Nonmonotonic grid')
        row={'state':state,'rows':len(rows)}
        for key in ['sc2c_pdos','zrcl2_pdos','total_projected_pdos']:
            row[key+'_at_EF_states_per_eV_cell']=interp(rows,key,0.0)
            row[key+'_window_weight_states_per_cell']=integrate(rows,key)
        out.append(row)
    return out


def read_crossings(root):
    # The .gnu energies are eV; use the parent SCF Fermi energy of this historical chain.
    text=(root/'scf/pwx.out').read_text()
    ef=float(re.findall(r'the Fermi energy is\s+([-+0-9.]+)',text)[-1])
    raw=[tuple(map(float,l.split())) for l in (root/'scf/bands.dat.gnu').read_text().splitlines() if l.strip()]
    lines=[l for l in (root/'scf/fatbands.projwfc_up').read_text().splitlines() if l.strip()]
    heads=[i for i,l in enumerate(lines[:30]) if len(l.split())==3 and all(x.isdigit() for x in l.split())]
    if len(heads)!=1: raise ValueError('Projection header not unique')
    idx=heads[0];nw,nk,nb=map(int,lines[idx].split())
    if (nw,nk,nb)!=(45,151,31) or lines[idx+1].split()!=['F','F']: raise ValueError('Wrong model/spin shape')
    if len(raw)!=nk*nb: raise ValueError('Band shape mismatch')
    bands=[raw[b*nk:(b+1)*nk] for b in range(nb)]
    if any([x for x,e in b] != [x for x,e in bands[0]] for b in bands): raise ValueError('Path mismatch')
    elements={1:'Zr',2:'C',3:'Cl',4:'Cl',5:'Sc',6:'Sc'}
    channels={'Zr_d':(1,'D'),'Sc_d':(5,'D'),'C_p':(2,'P'),'Cl_p':(3,'P')}
    w={k:[[0.0]*nk for _ in range(nb)] for k in [*channels,'all_projectors']}
    ptr=idx+2
    for expected in range(1,nw+1):
        h=lines[ptr].split();ptr+=1
        sid,atom=int(h[0]),int(h[1]);el=h[2];orb=h[3].upper()
        if sid!=expected or elements.get(atom)!=el: raise ValueError('State identity mismatch')
        angular=[c for c in orb if c in 'SPDF']
        if len(angular)!=1: raise ValueError('Invalid orbital label')
        c={'Zr':'Zr_d','Sc':'Sc_d','C':'C_p','Cl':'Cl_p'}[el]
        selected=angular[0]==('D' if el in ['Zr','Sc'] else 'P')
        for ik in range(nk):
            for ib in range(nb):
                a,b,v=lines[ptr].split();ptr+=1;v=float(v)
                if (int(a),int(b))!=(ik+1,ib+1) or not math.isfinite(v) or v < -1e-6:
                    raise ValueError('Projection row mismatch')
                w['all_projectors'][ib][ik]+=v
                if selected:w[c][ib][ik]+=v
    if ptr!=len(lines):raise ValueError('Trailing projection data')
    # Sum is diagnostic only. No complement is identified as an interstitial orbital.
    crossings=[]
    for ib,band in enumerate(bands):
        for ik,((xa,ea),(xb,eb)) in enumerate(zip(band,band[1:])):
            a,b=ea-ef,eb-ef
            if a*b < 0:
                f=-a/(b-a)
                row={'band':ib+1,'left_k_index':ik+1,'right_k_index':ik+2,
                     'path_distance_tpiba':xa+(xb-xa)*f,'EF_eV':ef}
                for c,arr in w.items():row[c]=arr[ib][ik]*(1-f)+arr[ib][ik+1]*f
                crossings.append(row)
    diagnostic={'shape':[nw,nk,nb],'EF_eV':ef,'all_projector_sum_min':min(v for b in w['all_projectors'] for v in b),
                'all_projector_sum_max':max(v for b in w['all_projectors'] for v in b),
                'crossing_method':'Linear interpolation of adjacent path energy and weights; not a uniform BZ integral.'}
    return crossings,diagnostic


def write_csv(p, rows):
    if not rows: raise ValueError('No result rows')
    with p.open('w',newline='') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)


def main():
    a=argparse.ArgumentParser();a.add_argument('data_root',type=Path);a.add_argument('--output',type=Path,default=Path('near-fermi-results'))
    args=a.parse_args();frozen=read_frozen(args.data_root);cross,checks=read_crossings(args.data_root)
    args.output.mkdir(parents=True,exist_ok=True)
    write_csv(args.output/'frozen-window.csv',frozen);write_csv(args.output/'path-crossings.csv',cross)
    report={'frozen_window_eV':[-0.1,0.0],'frozen':frozen,'historical_path':checks,'crossings':cross,
            'definitions':'Window integral is broadened projected spectral weight. It is not charge transfer or carrier density. Frozen QE7.2 and historical QE7.1 path chains are analyzed separately.'}
    (args.output/'near-fermi.json').write_text(json.dumps(report,indent=2)+'\n')
    print('Frozen cases:',len(frozen),'Path crossings:',len(cross))
    for r in frozen:
        print(r['state'], 'Sc2C D(EF)=',format(r['sc2c_pdos_at_EF_states_per_eV_cell'],'.6f'),
              'W[-0.1,0]=',format(r['sc2c_pdos_window_weight_states_per_cell'],'.6f'))
    for r in cross:print('Band',r['band'],'k bracket',r['left_k_index'],r['right_k_index'],
                         'Zr_d=',format(r['Zr_d'],'.6f'),'Sc_d=',format(r['Sc_d'],'.6f'))


if __name__=='__main__':main()

#!/usr/bin/env python3
from pathlib import Path
import csv
import re
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

ROOT = Path(__file__).resolve().parent.parent
OUT = ROOT / "delivery"
OUT.mkdir(exist_ok=True)
CASES = [
    ("hetero_0", "hetero_0", "Heterostructure 0%", 0.0, True),
    ("0", "sc2c_frozen_0", "Isolated Sc2C 0%", 0.0, False),
    ("hetero_0015", "hetero_0015", "Heterostructure +1.5%", 1.5, True),
    ("0015", "sc2c_frozen_0015", "Isolated Sc2C +1.5%", 1.5, False),
]
WINDOW = (-2.0, 2.0)
def load_case(spec):
    folder, prefix, name, strain, is_hetero = spec
    d = ROOT / folder
    pw = (d / "pwx.out").read_text(errors="replace")
    pr = (d / "projwfc.out").read_text(errors="replace")
    if "convergence has been achieved" not in pw:
        raise ValueError("SCF not converged: " + name)
    if "JOB DONE." not in pw or "JOB DONE." not in pr:
        raise ValueError("QE output incomplete: " + name)
    fm = next((x for x in pw.splitlines()
               if "the Fermi energy is" in x), None)
    if fm is None:
        raise ValueError("Fermi energy missing: " + name)
    ef = float(fm.split()[-2])
    total = np.loadtxt(d / (prefix + ".pdos_tot"), comments="#")
    energy = total[:, 0] - ef
    sc = np.zeros(len(energy))
    zr = np.zeros(len(energy))
    files = sorted(d.glob(prefix + ".pdos_atm#*_wfc#*"))
    expected = 19 if is_hetero else 10
    if len(files) != expected:
        raise ValueError(name + ": wrong projector file count")
    for f in files:
        match = re.search(r"atm#([0-9]+)", f.name)
        if match is None:
            raise ValueError("Bad projector filename: " + f.name)
        atom = int(match.group(1))
        part = np.loadtxt(f, comments="#")
        if len(part) != len(total):
            raise ValueError("PDOS row count differs: " + f.name)
        if not np.allclose(part[:, 0], total[:, 0], atol=1e-9):
            raise ValueError("PDOS energy grid differs: " + f.name)
        value = part[:, 2:].sum(axis=1)
        if not is_hetero or atom in (2, 5, 6):
            sc += value
        elif atom in (1, 3, 4):
            zr += value
        else:
            raise ValueError("Unexpected heterostructure atom")
    diff = sc + zr - total[:, 2]
    mask = np.abs(energy) <= 0.1
    area = np.trapz(np.abs(total[mask, 2]), energy[mask])
    l1 = 100 * np.trapz(np.abs(diff[mask]), energy[mask]) / area
    peak = np.max(np.abs(total[:, 2]))
    maxerr = 100 * np.max(np.abs(diff)) / peak
    return {"name": name, "strain": strain,
            "hetero": is_hetero, "ef": ef, "energy": energy,
            "sc": sc, "zr": zr, "total": total[:, 2],
            "rows": len(total), "step": np.median(np.diff(total[:, 0])),
            "nproj": len(files), "l1": l1, "maxerr": maxerr,
            "diff": diff}
states = [load_case(item) for item in CASES]
with (OUT / "frozen_pdos_long.csv").open("w", newline="") as f:
    writer = csv.writer(f)
    writer.writerow(["state", "strain_percent", "fermi_eV",
                     "energy_minus_fermi_eV", "sc2c_pdos",
                     "zrcl2_pdos", "total_projected_pdos",
                     "projected_sum_error"])
    for s in states:
        for i, x in enumerate(s["energy"]):
            writer.writerow([s["name"], s["strain"], s["ef"], x,
                             s["sc"][i], s["zr"][i],
                             s["total"][i], s["diff"][i]])
with (OUT / "frozen_pdos_summary.csv").open("w", newline="") as f:
    writer = csv.writer(f)
    writer.writerow(["state", "fermi_eV", "rows", "grid_step_eV",
                     "projector_files", "L1_sum_error_pct_EF_pm_0p1eV",
                     "max_sum_error_pct_total_peak"])
    for s in states:
        writer.writerow([s["name"], s["ef"], s["rows"], s["step"],
                         s["nproj"], s["l1"], s["maxerr"]])
fig, axes = plt.subplots(1, 2, figsize=(11, 6), sharex=True,
                         sharey=True)
blue, orange = "#0072B2", "#D55E00"
peaks = []
for strain, ax in zip((0.0, 1.5), axes):
    h = next(s for s in states if s["hetero"]
             and s["strain"] == strain)
    iso = next(s for s in states if not s["hetero"]
               and s["strain"] == strain)
    mask = (h["energy"] >= WINDOW[0]) & (h["energy"] <= WINDOW[1])
    mi = (iso["energy"] >= WINDOW[0]) & (iso["energy"] <= WINDOW[1])
    ax.plot(h["energy"][mask], h["sc"][mask], color=blue,
            lw=1.8, label="Sc2C layer: C#2 + Sc#5-6")
    ax.plot(iso["energy"][mi], iso["sc"][mi], color=blue,
            lw=1.8, ls="--", label="isolated Sc2C: Sc#1-2 + C#3")
    ax.plot(h["energy"][mask], h["zr"][mask], color=orange,
            lw=1.8, label="ZrCl2 layer: Zr#1 + Cl#3-4")
    peaks.extend([np.max(h["sc"][mask]),
                  np.max(iso["sc"][mi]), np.max(h["zr"][mask])])
    ax.axvline(0, color="#555555", lw=0.8, ls=":")
    ax.set_xlim(*WINDOW)
    ax.set_title("0% strain" if strain == 0 else "+1.5% strain")
    ax.set_xlabel("Energy relative to each Fermi level (eV)")
    ax.grid(axis="y", color="#DDDDDD", lw=0.6)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
for ax in axes:
    ax.set_ylim(0, max(peaks) * 1.05)
axes[0].set_ylabel("Projected DOS (states / eV / simulation cell)")
fig.suptitle("Frozen-cell projected DOS near the Fermi level",
             y=0.985, fontsize=14)
fig.text(0.5, 0.94,
         "Atom index = QE input order: hetero #1 Zr, #2 C, #3-4 Cl, #5-6 Sc.",
         ha="center", fontsize=9)
handles, labels = axes[0].get_legend_handles_labels()
fig.legend(handles, labels, loc="upper center",
           bbox_to_anchor=(0.5, 0.895), ncol=3,
           frameon=False, fontsize=8)
fig.text(0.5, 0.02,
         "Frozen coordinates, not relaxed equilibrium. EF alignment gives no absolute band offsets or charge transfer.",
         ha="center", fontsize=8, color="#444444")
fig.tight_layout(rect=[0, 0.07, 1, 0.83])
fig.savefig(OUT / "frozen_pdos.png", dpi=300, bbox_inches="tight")
fig.savefig(OUT / "frozen_pdos.pdf", bbox_inches="tight")
plt.close(fig)
print("Wrote products to", OUT)
for s in states:
    print(s["name"], "EF", s["ef"], "rows", s["rows"],
          "projectors", s["nproj"], "L1 error %", s["l1"],
          "max error %", s["maxerr"])