费米面几何嵌套

#!/usr/bin/env python3
"""Resolve original Al J(q) into band pairs without phonon matrix elements."""
from pathlib import Path
import argparse, csv, json
import numpy as np

def write_csv(path, rows):
    with path.open("w", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=list(rows[0]))
        writer.writeheader()
        writer.writerows(rows)

def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--al-root", type=Path, required=True,
                    help="Extracted al directory containing fermi/k24-cg and k32-cg")
    ap.add_argument("--output-dir", type=Path, required=True)
    args = ap.parse_args()
    if args.output_dir.exists():
        raise FileExistsError("Use a new output directory: " + str(args.output_dir))
    pairs, summaries, checks = [], [], []
    for n in [24, 32]:
        folder = args.al_root / "fermi" / f"k{n}-cg"
        with np.load(folder / "fermi-grid.npz", allow_pickle=False) as saved:
            energies = saved["energy_eV"].copy()
            grid = int(saved["grid"])
            ef = float(saved["fermi_eV"])
        info = json.loads((folder / "grid-info.json").read_text())
        if grid != n or energies.shape != (n, n, n, 6):
            raise ValueError("Expected the archived six-band complete Al grid")
        if not np.all(np.isfinite(energies)) or not np.isfinite(ef):
            raise ValueError("Non-finite energies or Fermi reference")
        if info["nks"] != n**3 or abs(info["fermi_eV"] - ef) > 1e-10:
            raise ValueError("Grid metadata differs")
        for sigma in [0.10, 0.20]:
            band_weight = np.exp(-0.5 * (energies / sigma)**2) / (
                sigma * np.sqrt(2 * np.pi))
            total_weight = band_weight.sum(axis=3)
            transform = np.fft.fftn(total_weight)
            joint = np.fft.ifftn(transform.conj() * transform).real / n**3
            zero = float(joint[0, 0, 0])
            tolerance = 1e-10 * max(1.0, zero)
            if zero <= 0 or np.min(joint) < -tolerance:
                raise ValueError("Invalid joint weight")
            if np.max(joint) > zero + tolerance:
                raise ValueError("Periodic autocorrelation exceeds J(0)")
            if abs(zero - float(np.mean(total_weight**2))) > tolerance:
                raise ValueError("q=0 direct sum differs")
            cut = np.loadtxt(folder / f"nesting-GX-s{sigma:.2f}.csv",
                             delimiter=",", skiprows=1)
            indices = np.arange(n//2 + 1)
            computed = joint[indices, 0, indices]
            if cut.shape != (n//2 + 1, 3):
                raise ValueError("Unexpected stored cut shape")
            if not np.allclose(cut[:, 0], indices/n, rtol=0, atol=1e-12):
                raise ValueError("Stored q coordinates differ")
            if not np.allclose(cut[:, 1], computed, rtol=1e-10, atol=1e-12):
                raise ValueError("Stored absolute J differs from FFT")
            if not np.allclose(cut[:, 2], computed/zero, rtol=1e-10, atol=1e-12):
                raise ValueError("Stored normalized J differs")
            for name, i in [("Gamma", 0), ("quarter", n//4), ("X", n//2)]:
                shift = (i, 0, i)
                shifted = np.roll(band_weight, tuple(-v for v in shift),
                                  axis=(0, 1, 2))
                matrix = (band_weight.reshape(-1, 6).T @
                          shifted.reshape(-1, 6)) / n**3
                direct = float(np.mean(total_weight *
                                       np.roll(total_weight, tuple(-v for v in shift),
                                               axis=(0, 1, 2))))
                total = float(matrix.sum())
                if abs(total - direct) > tolerance or abs(total - joint[shift]) > tolerance:
                    raise ValueError("Band-pair, direct and FFT sums differ")
                diagonal = float(np.trace(matrix))
                offdiagonal = float(matrix[~np.eye(6, dtype=bool)].sum())
                for band_n in range(6):
                    for band_m in range(6):
                        pairs.append({"kmesh": n, "sigma_eV": sigma,
                                      "q_name": name, "q_fraction_b1_plus_b3": i/n,
                                      "initial_band": band_n+1, "final_band": band_m+1,
                                      "J_pair_eV_minus2": float(matrix[band_n, band_m])})
                summaries.append({"kmesh": n, "sigma_eV": sigma,
                                  "q_name": name, "q_fraction_b1_plus_b3": i/n,
                                  "J_eV_minus2": total, "same_band_eV_minus2": diagonal,
                                  "different_band_eV_minus2": offdiagonal,
                                  "J_over_J0": total/zero})
                if name == "X":
                    print(f'{n}^3 sigma={sigma:.2f} eV: J(X)={total:.8f}; '
                          f'same-band={diagonal:.8f}; different-band={offdiagonal:.8f} eV^-2')
            checks.append({"kmesh": n, "sigma_eV": sigma,
                           "fermi_eV": ef, "all_six_bands_used": True,
                           "complete_grid_points": n**3, "stored_cut_reproduced": True,
                           "nonnegative_and_J_not_above_J0": True,
                           "band_pair_direct_FFT_agree": True})
    report = {"source_root": str(args.al_root), "numpy_version": np.__version__,
              "definition": "J_nm(q)=mean_k delta_sigma(E_n(k)-EF) delta_sigma(E_m(k+q)-EF)",
              "units": "eV^-2; all six bands; no added spin factor",
              "scope": "Geometric band-pair weights, not orbital/layer-resolved EPC or static susceptibility",
              "checks": checks, "summaries": summaries}
    args.output_dir.mkdir(parents=True)
    write_csv(args.output_dir / "band-pairs.csv", pairs)
    write_csv(args.output_dir / "band-pair-summary.csv", summaries)
    (args.output_dir / "band-pair-check.json").write_text(json.dumps(report, indent=2) + "\n")
    print(f'{len(summaries)} q-point summaries; {len(pairs)} band-pair records; '
          'pair/direct/FFT sums, stored cuts and J(q)<=J(0) checked.')
    print("No matrix elements, mode assignment, occupation denominator or DFT calculation added.")

if __name__ == "__main__":
    main()

from pathlib import Path
import numpy as np,xml.etree.ElementTree as E,json,hashlib,re
HARTREE_EV=27.211386245988
root=Path(__file__).resolve().parent
summary=[]
for n in [24,32]:
 d=root/f'k{n}-cg';out=(d/'al.nscf.out').read_text();err=(d/'al.nscf.err').read_text()
 assert out.count('JOB DONE.')==1 and not err and 'Error in routine' not in out
 assert 'not converged' not in out.lower()
 xml=d/'data-file-schema.xml';doc=E.parse(xml).getroot();o=doc.find('output');band=o.find('band_structure');nk=int(band.findtext('nks'));nb=int(band.findtext('nbnd'))
 assert nk==n**3
 b=np.array([np.fromstring(o.findtext('basis_set/reciprocal_lattice/'+tag),sep=' ') for tag in ['b1','b2','b3']])
 cell=np.array([np.fromstring(o.findtext('atomic_structure/cell/'+tag),sep=' ') for tag in ['a1','a2','a3']])*0.529177210903
 ef=float(band.findtext('fermi_energy'))*HARTREE_EV
 energies=np.empty((n,n,n,nb));seen=np.zeros((n,n,n),dtype=int)
 for point in band.findall('ks_energies'):
  cart=np.fromstring(point.findtext('k_point'),sep=' ');frac=cart@np.linalg.inv(b)
  scaled=frac*n;assert np.max(np.abs(scaled-np.rint(scaled)))<1e-7
  i=tuple(np.rint(scaled).astype(int)%n);seen[i]+=1
  energies[i]=np.fromstring(point.findtext('eigenvalues'),sep=' ')*HARTREE_EV-ef
 assert np.all(seen==1)
 np.savez_compressed(d/'fermi-grid.npz',energy_eV=energies,fermi_eV=ef,cell_angstrom=cell,grid=n)
 ranges=[{'band':j+1,'min_eV':float(energies[:,:,:,j].min()),'max_eV':float(energies[:,:,:,j].max())} for j in range(nb)]
 crossing=[r['band'] for r in ranges if r['min_eV']<0<r['max_eV']]
 record={'kmesh':n,'nks':nk,'fermi_eV':ef,'crossing_bands':crossing,'band_ranges':ranges,'source_xml_sha256':hashlib.sha256(xml.read_bytes()).hexdigest(),'nscf_out_sha256':hashlib.sha256((d/'al.nscf.out').read_bytes()).hexdigest(),'definition':'Energy grid includes all k, no interpolation, E-EF in eV.'}
 (d/'grid-info.json').write_text(json.dumps(record,indent=2));summary.append(record)
 print(f'k={n}^3 nks={nk} EF={ef:.8f} eV crossing bands={crossing}; all grid cells assigned once')
 for sigma in [.10,.20]:
  weight=np.exp(-0.5*(energies/sigma)**2).sum(axis=3)/(sigma*np.sqrt(2*np.pi))
  spectrum=np.fft.fftn(weight);J=np.fft.ifftn(spectrum.conj()*spectrum).real/nk
  assert np.min(J)>-1e-10
  assert abs(J[0,0,0]-np.mean(weight*weight))<1e-8
  # one nontrivial point cross-check against direct Brillouin-zone sum
  idx=(n//4,0,n//4);direct=np.mean(weight*np.roll(weight,tuple(-x for x in idx),axis=(0,1,2)))
  assert abs(J[idx]-direct)<1e-8
  np.savez_compressed(d/f'nesting-s{sigma:.2f}.npz',nesting_eV_minus2=J,sigma_eV=sigma,grid=n)
  cut=np.array([[i/n,J[i,0,i],J[i,0,i]/J[0,0,0]] for i in range(n//2+1)])
  np.savetxt(d/f'nesting-GX-s{sigma:.2f}.csv',cut,delimiter=',',header='q_fraction_along_b1_plus_b3,J_eV_minus2,J_over_J0',comments='')
  print(f' sigma={sigma:.2f} eV J(0)={J[0,0,0]:.8f} J(X)={J[n//2,0,n//2]:.8f} eV^-2; direct-sum check passed')
(root/'summary.json').write_text(json.dumps(summary,indent=2))

from atlas_plot_style import install as install_atlas_style
install_atlas_style()
from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt
r=Path(__file__).resolve().parent
fig,axes=plt.subplots(1,2,figsize=(10.5,4.2),layout="constrained")
for n in [24,32]:
    for sigma in [.10,.20]:
        x=np.loadtxt(r/f"fermi/k{n}-cg/nesting-GX-s{sigma:.2f}.csv",delimiter=",",skiprows=1)
        label=f"{n}³, σ={sigma:.2f} eV"
        for ax,col in zip(axes,[1,2]):ax.plot(x[:,0],x[:,col],"o-",ms=3,lw=1.3,label=label)
for ax in axes:ax.set(xlabel="q = t(b₁+b₃), Γ → X",xlim=(0,.5));ax.grid(alpha=.2)
axes[0].set_ylabel("J(q) (eV⁻²)");axes[1].set_ylabel("J(q) / J(0)");axes[1].legend(frameon=False,fontsize=8)
(r/"figures").mkdir(exist_ok=True)
fig.savefig(r/"figures/fermi-nesting.png",dpi=220);fig.savefig(r/"figures/fermi-nesting.pdf")

# Run from the extracted al/ directory: gnuplot plot_nesting.gnu
if (!exists("data_root")) data_root="."
if (!exists("output_path")) output_path="fermi-nesting-gnuplot.png"
set encoding utf8
set datafile separator comma
file(n,s)=sprintf("%s/fermi/k%d-cg/nesting-GX-s%.2f.csv",data_root,n,s)
set terminal pngcairo enhanced font "DejaVu Sans,12" size 1500,600
set output output_path
set multiplot layout 1,2 margins 0.075,0.98,0.17,0.83 spacing 0.10
set xrange [0:0.5]
set yrange [0:*]
set xlabel "q = t(b_1+b_3), Gamma to X"
set grid ytics lc rgb "#dddddd"
set tics nomirror
set key top right font ",10"
set title "(a) Absolute joint weight"
set ylabel "J(q) (eV^{-2})"
plot file(24,0.10) using 1:2 with linespoints lw 1.5 pt 7 ps 0.55 lc rgb "#0072b2" title "24^3, sigma=0.10 eV", \
 file(24,0.20) using 1:2 with linespoints lw 1.5 pt 5 ps 0.55 lc rgb "#d55e00" title "24^3, sigma=0.20 eV", \
 file(32,0.10) using 1:2 with linespoints lw 1.5 pt 9 ps 0.55 lc rgb "#009e73" title "32^3, sigma=0.10 eV", \
 file(32,0.20) using 1:2 with linespoints lw 1.5 pt 11 ps 0.55 lc rgb "#cc79a7" title "32^3, sigma=0.20 eV"
set title "(b) Shape normalized to q=0"
set ylabel "J(q) / J(0)"
plot file(24,0.10) using 1:3 with linespoints lw 1.5 pt 7 ps 0.55 lc rgb "#0072b2" title "24^3, sigma=0.10 eV", \
 file(24,0.20) using 1:3 with linespoints lw 1.5 pt 5 ps 0.55 lc rgb "#d55e00" title "24^3, sigma=0.20 eV", \
 file(32,0.10) using 1:3 with linespoints lw 1.5 pt 9 ps 0.55 lc rgb "#009e73" title "32^3, sigma=0.10 eV", \
 file(32,0.20) using 1:3 with linespoints lw 1.5 pt 11 ps 0.55 lc rgb "#cc79a7" title "32^3, sigma=0.20 eV"
unset multiplot
print "Read the four original CSVs; no interpolation, smoothing or peak normalization beyond stored column3 J/J(0)."