#!/usr/bin/env amspython from __future__ import annotations from dataclasses import dataclass from pathlib import Path import shutil import matplotlib.pyplot as plt import pandas as pd from scm.base import Units from scm.plams import AMSJob ROOT = Path(__file__).resolve().parent WORKDIR = ROOT / "01-run_workdir" REPORT_MD = ROOT / "report.md" REPORT_BK = ROOT / "report.md.bk" PLOT_PNG = ROOT / "relative_energy_vs_interlayer_distance.png" N_ATOMS = 4 HARTREE_TO_MEV = Units.conversion_factor("hartree", "eV") * 1000.0 @dataclass(frozen=True) class JobRef: name: str label: str justification: str JOB_REFS: tuple[JobRef, ...] = ( JobRef("pbe_psl_paw_seg1", "PBE / pslibrary-PAW, 6.0-9.0 A", "Short-range scan for PBE with pslibrary-PAW."), JobRef("pbe_psl_paw_seg2", "PBE / pslibrary-PAW, 9.0-16.0 A", "Long-range continuation for PBE with pslibrary-PAW."), JobRef( "b86bpbe_xdm_psl_paw_seg1", "B86bPBE-XDM / pslibrary-PAW, 6.0-9.0 A", "Short-range scan for B86bPBE with QE XDM dispersion and pslibrary-PAW.", ), JobRef( "b86bpbe_xdm_psl_paw_seg2", "B86bPBE-XDM / pslibrary-PAW, 9.0-16.0 A", "Long-range continuation for B86bPBE with QE XDM dispersion and pslibrary-PAW.", ), JobRef("r2scan_dojo_seg1", "r2SCAN / Dojo, 6.0-9.0 A", "Short-range scan for r2SCAN with Dojo pseudopotentials."), JobRef("r2scan_dojo_seg2", "r2SCAN / Dojo, 9.0-16.0 A", "Long-range continuation for r2SCAN with Dojo pseudopotentials."), JobRef( "r2scan_d4_dojo_seg1", "r2SCAN-D4 / Dojo, 6.0-9.0 A", "Short-range scan for r2SCAN with Dojo pseudopotentials and the AMS D4 add-on.", ), JobRef( "r2scan_d4_dojo_seg2", "r2SCAN-D4 / Dojo, 9.0-16.0 A", "Long-range continuation for r2SCAN with Dojo pseudopotentials and the AMS D4 add-on.", ), ) def load_job(job_name: str) -> AMSJob: return AMSJob.load_external(WORKDIR / job_name) def build_curve(job_names: tuple[str, str], label: str) -> pd.DataFrame: frames: list[pd.DataFrame] = [] for job_name in job_names: job = load_job(job_name) pes = job.results.get_pesscan_results(molecules=False) c_index = pes["RaveledScanCoords"].index("c") c_unit = pes["RaveledUnits"][c_index] c_conversion = Units.conversion_factor(c_unit, "angstrom") c_values = pd.Series(pes["RaveledPESCoords"][c_index], dtype=float) * c_conversion energies = pd.Series(pes["PES"], dtype=float) frames.append(pd.DataFrame({"c_angstrom": c_values, "energy_hartree": energies})) curve = pd.concat(frames, ignore_index=True) curve = curve.drop_duplicates(subset="c_angstrom", keep="first").sort_values("c_angstrom").reset_index(drop=True) reference_energy = curve.loc[curve["c_angstrom"].idxmax(), "energy_hartree"] curve["c_over_2_angstrom"] = curve["c_angstrom"] / 2.0 curve["relative_mev_per_atom"] = (curve["energy_hartree"] - reference_energy) * HARTREE_TO_MEV / N_ATOMS curve["functional"] = label curve["minimum_relative_mev_per_atom"] = curve["relative_mev_per_atom"].min() return curve def plot_curves(curves: list[pd.DataFrame]) -> None: curves_to_plot = sorted(curves, key=lambda df: float(df["minimum_relative_mev_per_atom"].iloc[0]), reverse=True) plt.figure(figsize=(7.2, 4.8)) for curve in curves_to_plot: label = str(curve["functional"].iloc[0]) plt.plot( curve["c_over_2_angstrom"], curve["relative_mev_per_atom"], marker="o", markersize=3, linewidth=1.4, label=label, ) plt.xlabel("c/2 (angstrom)") plt.ylabel("Relative energy (meV/atom)") plt.legend(frameon=False) plt.tight_layout() plt.savefig(PLOT_PNG, dpi=200) plt.close() def build_summary_table(curves: list[pd.DataFrame]) -> pd.DataFrame: rows: list[dict[str, float | str]] = [] for curve in curves: minimum_row = curve.loc[curve["relative_mev_per_atom"].idxmin()] rows.append( { "Functional": str(curve["functional"].iloc[0]), "Min at c/2 (A)": float(minimum_row["c_over_2_angstrom"]), "Min rel. energy (meV/atom)": float(minimum_row["relative_mev_per_atom"]), "Ref. point c/2 (A)": float(curve["c_over_2_angstrom"].max()), } ) table = pd.DataFrame(rows) return table.sort_values("Min rel. energy (meV/atom)").reset_index(drop=True) def render_report(curves: list[pd.DataFrame]) -> None: if REPORT_MD.exists(): shutil.copy2(REPORT_MD, REPORT_BK) summary_table = build_summary_table(curves) provenance_lines: list[str] = [] input_sections: list[str] = [] for job_ref in JOB_REFS: job = load_job(job_ref.name) provenance_lines.append(f"- `{job_ref.name}`: {job_ref.justification}") input_sections.append(f"### `{job_ref.name}`\n\nReason: {job_ref.justification}\n\n```text\n{job.get_input().strip()}\n```") report_text = "\n".join( [ "# Graphite c-axis PES scans with QuantumESPRESSO", "", "This report combines two non-optimizing PES scan segments per functional for graphite, scanning the c lattice parameter and plotting the relative energy against the interlayer distance c/2.", "", "## Conclusion", "", "All curves are referenced to the energy at the longest sampled c value for each functional. The plot order is the reverse order of minimum energy so the legend follows the visible stacking order of the plotted curves.", "", "## Relative-energy plot", "", f"![Relative energy versus interlayer distance]({PLOT_PNG.name})", "", "## Summary table", "", summary_table.to_markdown(index=False, floatfmt=".6f"), "", "## Provenance", "", *provenance_lines, "", "## Calculation inputs", "", *input_sections, "", ] ) REPORT_MD.write_text(report_text) def main() -> None: curves = [ build_curve(("pbe_psl_paw_seg1", "pbe_psl_paw_seg2"), "PBE / pslibrary-PAW"), build_curve(("b86bpbe_xdm_psl_paw_seg1", "b86bpbe_xdm_psl_paw_seg2"), "B86bPBE-XDM / pslibrary-PAW"), build_curve(("r2scan_dojo_seg1", "r2scan_dojo_seg2"), "r2SCAN / Dojo"), build_curve(("r2scan_d4_dojo_seg1", "r2scan_d4_dojo_seg2"), "r2SCAN-D4 / Dojo"), ] plot_curves(curves) render_report(curves) if __name__ == "__main__": main()