应变下的电荷、电子态与声子/EPC

# Read the published CSV directly; all energies are already relative to each EF.
# Run: gnuplot plot_frozen_pdos.gnu
# Optional: gnuplot -e "energy_window=0.5" plot_frozen_pdos.gnu
if (!exists("datafile")) datafile="frozen_pdos_long.csv"
if (!exists("output_path")) output_path="frozen-pdos-gnuplot.png"
if (!exists("energy_window")) energy_window=2.0
set encoding utf8
set datafile separator comma
# Shared y range from every plotted Sc2C/ZrCl2 trace within the chosen window.
stats datafile using (abs($4)<=energy_window ? $5 : 1/0) nooutput
shared_max=STATS_max
stats datafile using (abs($4)<=energy_window ? $6 : 1/0) nooutput
shared_max=1.08*(shared_max>STATS_max ? shared_max : STATS_max)
set terminal pngcairo enhanced font "DejaVu Sans,12" size 1500,620
set output output_path
set multiplot layout 1,2 margins 0.075,0.97,0.18,0.78 spacing 0.06
set xrange [-energy_window:energy_window]
set yrange [0:shared_max]
set xlabel "Energy relative to each Fermi level (eV)"
set ylabel "Projected DOS (states / eV / simulation cell)"
set grid ytics lc rgb "#dddddd"
set border 3
set tics nomirror
set key top right font ",10"
set arrow 1 from 0,graph 0 to 0,graph 1 nohead dt 2 lc rgb "#777777"
set title "(a) 0% strain | frozen geometry"
plot datafile using 4:(stringcolumn(1) eq "Heterostructure 0%" ? $5 : 1/0) with lines lw 2 lc rgb "#0072b2" title "Sc2C layer", \
 datafile using 4:(stringcolumn(1) eq "Isolated Sc2C 0%" ? $5 : 1/0) with lines lw 2 dt 2 lc rgb "#0072b2" title "matched isolated Sc2C", \
 datafile using 4:(stringcolumn(1) eq "Heterostructure 0%" ? $6 : 1/0) with lines lw 2 lc rgb "#d55e00" title "ZrCl2 layer"
set title "(b) +1.5% strain | frozen geometry"
unset ylabel
plot datafile using 4:(stringcolumn(1) eq "Heterostructure +1.5%" ? $5 : 1/0) with lines lw 2 lc rgb "#0072b2" title "Sc2C layer", \
 datafile using 4:(stringcolumn(1) eq "Isolated Sc2C +1.5%" ? $5 : 1/0) with lines lw 2 dt 2 lc rgb "#0072b2" title "matched isolated Sc2C", \
 datafile using 4:(stringcolumn(1) eq "Heterostructure +1.5%" ? $6 : 1/0) with lines lw 2 lc rgb "#d55e00" title "ZrCl2 layer"
unset multiplot
print sprintf("Read column4 and layer columns5/6 without resampling; |E-EF|<=%.2f eV; shared_ymax=%.8f", energy_window,shared_max)

#!/usr/bin/env python3
"""Inspect accepted frozen-control PDOS tables; no DFT or plotting."""
from pathlib import Path
import csv, json, math, sys
root = Path(__file__).resolve().parent
source = Path(sys.argv[1]) if len(sys.argv) > 1 else root / "frozen_pdos_long.csv"
expected = ["Heterostructure 0%", "Isolated Sc2C 0%",
            "Heterostructure +1.5%", "Isolated Sc2C +1.5%"]
groups = {name: [] for name in expected}
with source.open(newline="") as handle:
    reader = csv.DictReader(handle)
    names = ["strain_percent", "fermi_eV", "energy_minus_fermi_eV",
             "sc2c_pdos", "zrcl2_pdos", "total_projected_pdos",
             "projected_sum_error"]
    if not set(["state"] + names).issubset(reader.fieldnames or []):
        raise ValueError("Missing PDOS fields")
    for item in reader:
        state = item["state"]
        if state not in groups:
            raise ValueError("Unknown state: " + state)
        values = {k: float(item[k]) for k in names}
        if not all(math.isfinite(v) for v in values.values()):
            raise ValueError("Non-finite PDOS row")
        discrepancy = values["sc2c_pdos"] + values["zrcl2_pdos"] - values["total_projected_pdos"]
        if abs(discrepancy - values["projected_sum_error"]) > 1e-10:
            raise ValueError("Saved projection mismatch differs")
        groups[state].append(values)
rows, checks = [], []
for state in expected:
    data = groups[state]
    if not data:
        raise ValueError("Empty state: " + state)
    if len({x["strain_percent"] for x in data}) != 1 or len({x["fermi_eV"] for x in data}) != 1:
        raise ValueError("Mixed strain or Fermi reference")
    energy = [x["energy_minus_fermi_eV"] for x in data]
    steps = [b-a for a,b in zip(energy[:-1],energy[1:])]
    if not steps or min(steps) <= 0 or max(abs(x-0.005) for x in steps) > 1e-8:
        raise ValueError("Unexpected energy grid")
    if not min(energy) < 0 < max(energy):
        raise ValueError("Grid does not bracket EF")
    near = min(data, key=lambda x: abs(x["energy_minus_fermi_eV"]))
    if abs(near["energy_minus_fermi_eV"]) > 0.0025 + 1e-8:
        raise ValueError("No near-Fermi grid point")
    window = [x for x in data if abs(x["energy_minus_fermi_eV"]) <= 0.1]
    denominator = sum(abs(x["total_projected_pdos"]) for x in window)
    if not window or denominator <= 0:
        raise ValueError("Invalid closure window")
    closure = 100 * sum(abs(x["projected_sum_error"]) for x in window) / denominator
    row = {"state": state, "strain_percent": near["strain_percent"],
           "nearest_E_minus_EF_eV": near["energy_minus_fermi_eV"],
           "sc2c_pdos": near["sc2c_pdos"], "zrcl2_pdos": near["zrcl2_pdos"],
           "total_projected_pdos": near["total_projected_pdos"]}
    rows.append(row)
    checks.append({"state":state, "rows":len(data), "grid_step_eV":steps[0],
                   "L1_projection_error_percent_EF_pm_0p1":closure})
    print(f'{state}: E-EF={row["nearest_E_minus_EF_eV"]:+.4f} eV; '
          f'Sc2C={row["sc2c_pdos"]:.6f}; ZrCl2={row["zrcl2_pdos"]:.6f}; '
          f'total={row["total_projected_pdos"]:.2f} states/(eV cell)')
with (root/"nearest-fermi-pdos.csv").open("w", newline="") as handle:
    writer = csv.DictWriter(handle, fieldnames=list(rows[0]))
    writer.writeheader(); writer.writerows(rows)
result = {"source":source.name, "method":"Nearest grid point; no interpolation at EF",
          "units":"PDOS: states/(eV cell); E-EF: eV",
          "reference":"Each state has its own Fermi energy",
          "scope":"Frozen geometry; projected DOS is not transferred charge",
          "rows":rows, "checks":checks}
(root/"nearest-fermi-pdos.json").write_text(json.dumps(result,indent=2)+"\n")
print("Four frozen-control PDOS states checked; no DFT invoked.")

#!/usr/bin/env python3
"""Integrate a fixed own-EF PDOS window from the published frozen controls."""
from pathlib import Path
import argparse, bisect, csv, json, math

STATES = {
    "Heterostructure 0%": 0.0,
    "Isolated Sc2C 0%": 0.0,
    "Heterostructure +1.5%": 1.5,
    "Isolated Sc2C +1.5%": 1.5,
}
COLUMNS = ["sc2c_pdos", "zrcl2_pdos", "total_projected_pdos"]

def integrate(x, y, lo, hi):
    if not x[0] < lo < hi < x[-1]:
        raise ValueError("Window must lie inside the saved energy grid")
    def at(t):
        i = bisect.bisect_right(x, t) - 1
        return y[i] + (y[i+1] - y[i]) * (t - x[i]) / (x[i+1] - x[i])
    points = [(lo, at(lo))]
    points.extend((a, b) for a, b in zip(x, y) if lo < a < hi)
    points.append((hi, at(hi)))
    return math.fsum((b[0] - a[0]) * (a[1] + b[1]) / 2
                     for a, b in zip(points, points[1:]))

def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--input", type=Path, required=True)
    ap.add_argument("--half-width", type=float, default=0.1)
    ap.add_argument("--output-dir", type=Path, required=True)
    args = ap.parse_args()
    if not math.isfinite(args.half_width) or args.half_width <= 0:
        raise ValueError("half-width must be a positive finite eV value")
    if args.output_dir.exists():
        raise FileExistsError("Use a new output directory: " + str(args.output_dir))
    groups = {state: [] for state in STATES}
    fields = ["strain_percent", "fermi_eV", "energy_minus_fermi_eV",
              *COLUMNS, "projected_sum_error"]
    with args.input.open(newline="") as handle:
        reader = csv.DictReader(handle)
        if not {"state", *fields}.issubset(reader.fieldnames or []):
            raise ValueError("Missing PDOS fields")
        for item in reader:
            state = item["state"]
            if state not in groups:
                raise ValueError("Unknown state: " + state)
            row = {name: float(item[name]) for name in fields}
            if not all(math.isfinite(v) for v in row.values()):
                raise ValueError("Non-finite data: " + state)
            mismatch = row["sc2c_pdos"] + row["zrcl2_pdos"] - row["total_projected_pdos"]
            if abs(mismatch - row["projected_sum_error"]) > 1e-10:
                raise ValueError("Projection discrepancy column does not match")
            groups[state].append(row)
    rows, checks = [], []
    width = 2 * args.half_width
    for state, strain in STATES.items():
        data = groups[state]
        if not data or {v["strain_percent"] for v in data} != {strain}:
            raise ValueError("Empty or mixed-strain state: " + state)
        if len({v["fermi_eV"] for v in data}) != 1:
            raise ValueError("Mixed Fermi energies: " + state)
        x = [v["energy_minus_fermi_eV"] for v in data]
        steps = [b-a for a, b in zip(x, x[1:])]
        if not steps or any(abs(v - 0.005) > 1e-8 for v in steps):
            raise ValueError("Expected a strictly increasing 0.005 eV grid")
        weights = {col: integrate(x, [v[col] for v in data],
                                  -args.half_width, args.half_width)
                   for col in COLUMNS}
        absolute_mismatch = integrate(x, [abs(v["projected_sum_error"]) for v in data],
                                      -args.half_width, args.half_width)
        if weights["total_projected_pdos"] <= 0:
            raise ValueError("Nonpositive integrated total projection")
        row = {"state": state, "strain_percent": strain,
               "half_width_eV": args.half_width, "fermi_eV": data[0]["fermi_eV"]}
        for col in COLUMNS:
            row[col + "_weight_states_cell"] = weights[col]
            row[col + "_mean_states_eV_cell"] = weights[col] / width
        rows.append(row)
        checks.append({"state": state, "source_rows": len(data),
                       "L1_closure_percent": 100 * absolute_mismatch /
                       weights["total_projected_pdos"]})
        print(f'{state}: Sc2C={weights["sc2c_pdos"]:.6f}; '
              f'ZrCl2={weights["zrcl2_pdos"]:.6f}; '
              f'total projection={weights["total_projected_pdos"]:.6f} states/cell')
    changes = {}
    for tag, first, last in [("heterostructure", rows[0], rows[2]),
                              ("isolated_Sc2C", rows[1], rows[3])]:
        key = "sc2c_pdos_weight_states_cell"
        changes[tag] = 100 * (last[key] / first[key] - 1)
    print(f'Sc2C window-weight change: heterostructure={changes["heterostructure"]:+.4f}%; '
          f'isolated={changes["isolated_Sc2C"]:+.4f}%')
    report = {"source": str(args.input), "half_width_eV": args.half_width,
              "method": "Piecewise-linear spectrum; exact window endpoints; trapezoidal integral",
              "reference": "Each state uses its own saved Fermi energy; no absolute band alignment",
              "units": "Integrated projected spectral weight: states/cell; mean PDOS: states/(eV cell)",
              "scope": "Frozen geometry; no occupation integral, transferred charge, or EPC",
              "rows": rows, "checks": checks, "sc2c_weight_change_percent": changes}
    args.output_dir.mkdir(parents=True)
    with (args.output_dir / "window-pdos.csv").open("w", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=list(rows[0]))
        writer.writeheader()
        writer.writerows(rows)
    (args.output_dir / "window-pdos.json").write_text(json.dumps(report, indent=2) + "\n")
    print(f'Window [-{args.half_width:g}, +{args.half_width:g}] eV; '
          'linear endpoints only; no DFT or charge-transfer inference.')

if __name__ == "__main__":
    main()

# gnuplot 6.0; direct published CSV, piecewise-linear sample connections.
# Run beside frozen_pdos_long.csv: gnuplot plot_window_pdos.gnu
if (!exists("datafile")) datafile="frozen_pdos_long.csv"
if (!exists("output_base")) output_base="pdos-window"
set encoding utf8
set datafile separator comma
stats datafile using (abs($4)<=0.5 ? $5 : 1/0) nooutput
shared_max=STATS_max
stats datafile using (abs($4)<=0.5 ? $6 : 1/0) nooutput
shared_max=1.08*(shared_max>STATS_max ? shared_max : STATS_max)
do for [format_index=1:2] {
    if (format_index==1) {
        set terminal pngcairo enhanced font "DejaVu Sans,12" size 1500,620
        set output output_base.".png"
    } else {
        set terminal svg enhanced font "DejaVu Sans,12" size 1500,620
        set output output_base.".svg"
    }
    set multiplot layout 1,2 margins 0.075,0.97,0.18,0.78 spacing 0.06
    set xrange [-0.5:0.5]
    set yrange [0:shared_max]
    set xlabel "Energy relative to each Fermi level (eV)"
    set ylabel "Projected DOS (states / eV / simulation cell)"
    set grid ytics lc rgb "#dddddd"
    set border 3
    set tics nomirror
    set key top right font ",10"
    set object 1 rect from -0.1,graph 0 to 0.1,graph 1 behind \
        fc rgb "#777777" fs transparent solid 0.10 noborder
    set arrow 1 from 0,graph 0 to 0,graph 1 nohead dt 2 lc rgb "#777777"
    set arrow 2 from -0.1,graph 0 to -0.1,graph 1 nohead dt 3 lc rgb "#aaaaaa"
    set arrow 3 from 0.1,graph 0 to 0.1,graph 1 nohead dt 3 lc rgb "#aaaaaa"
    set label 1 "integrated window" at 0,graph 0.95 center font ",10" tc rgb "#555555"
    set title "(a) 0% strain | frozen geometry"
    plot datafile using 4:(stringcolumn(1) eq "Heterostructure 0%" ? $5 : 1/0) \
        with lines lw 2 lc rgb "#0072b2" title "Sc2C layer", \
      datafile using 4:(stringcolumn(1) eq "Isolated Sc2C 0%" ? $5 : 1/0) \
        with lines lw 2 dt 2 lc rgb "#0072b2" title "matched isolated Sc2C", \
      datafile using 4:(stringcolumn(1) eq "Heterostructure 0%" ? $6 : 1/0) \
        with lines lw 2 lc rgb "#d55e00" title "ZrCl2 layer"
    set title "(b) +1.5% strain | frozen geometry"
    unset ylabel
    plot datafile using 4:(stringcolumn(1) eq "Heterostructure +1.5%" ? $5 : 1/0) \
        with lines lw 2 lc rgb "#0072b2" title "Sc2C layer", \
      datafile using 4:(stringcolumn(1) eq "Isolated Sc2C +1.5%" ? $5 : 1/0) \
        with lines lw 2 dt 2 lc rgb "#0072b2" title "matched isolated Sc2C", \
      datafile using 4:(stringcolumn(1) eq "Heterostructure +1.5%" ? $6 : 1/0) \
        with lines lw 2 lc rgb "#d55e00" title "ZrCl2 layer"
    unset multiplot
    unset output
}
print sprintf("Original CSV, common y range 0..%.8f; window +/-0.1 eV; PNG/SVG written.",shared_max)

from pathlib import Path
import csv, re, json, hashlib, math
r=Path(__file__).resolve().parent
rows=[]; hashes={}
def cell(text):
    rows=re.findall(r'CELL_PARAMETERS\s+angstrom\s*\n([^\n]+)\n([^\n]+)\n([^\n]+)',text)
    if len(rows)!=1:raise ValueError('Expected one angstrom cell')
    return [[float(v) for v in line.split()] for line in rows[0]]
reference=cell((r/'elastic/al.reference.in').read_text())

for case in json.loads((r/'elastic/cases.json').read_text()):
    if case['mode']!='xx':continue
    p=r/'elastic'/case['label'];s=(p/'al.scf.out').read_text();err=(p/'al.scf.err').read_text()
    if s.count('JOB DONE.')!=1 or 'convergence has been achieved' not in s or err.strip() or 'convergence NOT achieved' in s or 'Error in routine' in s:
        raise ValueError('Incomplete SCF: '+case['label'])
    inp=(p/'al.scf.in').read_text(); current=cell(inp)
    strain=float(case['engineering_strain'])
    for i in range(3):
        for j in range(3):
            expected=reference[i][j]*(1+strain if j==0 else 1)
            if abs(current[i][j]-expected)>1e-10:raise ValueError('Cell/strain mismatch: '+case['label'])
    energy=float(re.findall(r'!\s+total energy\s+=\s+([-0-9.]+)',s)[-1])
    block=re.findall(r'total\s+stress[^\n]*\n([^\n]+)\n([^\n]+)\n([^\n]+)',s)[-1]
    stress=[[-14710.5076*float(x) for x in line.split()[:3]] for line in block]
    if not math.isfinite(energy) or not all(math.isfinite(v) for row in stress for v in row):raise ValueError('Non-finite result')
    rows.append({**case,'energy_Ry':energy,'sigma_xx_GPa':stress[0][0],'sigma_yy_GPa':stress[1][1],'sigma_zz_GPa':stress[2][2],'sigma_xy_GPa':stress[0][1]})
    for n in ['al.scf.in','al.scf.out','al.scf.err']:
        hashes[str((p/n).relative_to(r))]=hashlib.sha256((p/n).read_bytes()).hexdigest()
rows.sort(key=lambda x:x['engineering_strain'])
if len(rows)!=6:raise ValueError('Expected six xx samples')
with (r/'elastic/strain-stress.csv').open('w') as f:
    w=csv.DictWriter(f,fieldnames=rows[0]);w.writeheader();w.writerows(rows)
(r/'strain-summary.json').write_text(json.dumps({'rows':rows,'source_sha256':hashes,'energy':'QE printed F in Ry/one-atom cell','stress':'Tensile-positive; minus QE Ry/bohr^3 * 14710.5076 GPa','scope':'Six xx fixed-cell SCFs, 16^3 k mesh; no doping/DOS/EPC.'},indent=2)+'\n')
print('Six xx SCFs: electronic convergence, one JOB DONE, empty stderr.')
for x in rows:print(f"{x['engineering_strain']:+.3f} F={x['energy_Ry']:.8f} Ry sigma_xx={x['sigma_xx_GPa']:+.6f} sigma_yy={x['sigma_yy_GPa']:+.6f} GPa")

from atlas_plot_style import install as install_atlas_style
install_atlas_style()
from pathlib import Path
import csv
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
r=Path(__file__).resolve().parent
rows=sorted((x for x in csv.DictReader((r/'elastic/strain-stress.csv').open()) if x['mode']=='xx'),key=lambda x:float(x['engineering_strain']))
e=np.array([float(x['engineering_strain']) for x in rows])*100
energy=np.array([float(x['energy_Ry']) for x in rows])
sxx=[float(x['sigma_xx_GPa']) for x in rows]
syy=[float(x['sigma_yy_GPa']) for x in rows]
fig,axes=plt.subplots(1,2,figsize=(9,3.8),layout='constrained')
axes[0].plot(e,(energy-energy.min())*13.605693122994*1000,'o-',color='#0072b2')
axes[0].set(xlabel='Applied x strain (%)',ylabel='F - lowest sampled F (meV/atom)')
axes[1].plot(e,sxx,'o-',label='Tensile-positive stress xx',color='#0072b2')
axes[1].plot(e,syy,'s-',label='Tensile-positive stress yy',color='#d55e00')
axes[1].set(xlabel='Applied x strain (%)',ylabel='Stress (GPa)')
axes[1].legend(frameon=False,fontsize=8)
for ax in axes:
    ax.axvline(0,color='.6',lw=.7);ax.grid(alpha=.18)
fig.suptitle('fcc Al | fixed electron number, 16³ k mesh, MV width 0.02 Ry',fontsize=11)
(r/'figures').mkdir(exist_ok=True)
fig.savefig(r/'figures/strain-scan.png',dpi=220)
fig.savefig(r/'figures/strain-scan.svg')