#!/usr/bin/env amspython from __future__ import annotations import os import shutil import subprocess from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from scm.plams import AMSJob, view ROOT = Path(__file__).resolve().parent FIGURES = ROOT / "figures" TABLES = ROOT / "tables" def latest_workdir() -> Path: candidates = [path for path in ROOT.glob("01-run_workdir*") if path.is_dir()] if not candidates: raise FileNotFoundError("No 01-run_workdir directory was found") return max(candidates, key=lambda path: path.stat().st_mtime) def frontier_rows(job: AMSJob, method: str, count: int = 3) -> list[dict[str, object]]: energies = np.asarray(job.results.get_orbital_energies(unit="eV"))[0] occupations = np.asarray(job.results.get_orbital_occupations())[0] occupied = np.flatnonzero(occupations > 1.0e-6) virtual = np.flatnonzero(occupations <= 1.0e-6) if len(occupied) < count or len(virtual) < count: raise RuntimeError(f"Too few orbitals were returned for {method}") homo_index = int(occupied[-1]) lumo_index = int(virtual[0]) selected = list(occupied[-count:]) + list(virtual[:count]) rows: list[dict[str, object]] = [] for index in selected: if index <= homo_index: offset = homo_index - int(index) label = "HOMO" if offset == 0 else f"HOMO-{offset}" else: offset = int(index) - lumo_index label = "LUMO" if offset == 0 else f"LUMO+{offset}" rows.append( { "Method": method, "Orbital": label, "Energy (eV)": float(energies[index]), "Occupation": float(occupations[index]), } ) return rows def make_energy_level_plot(orbital_table: pd.DataFrame) -> None: fig, ax = plt.subplots(figsize=(5.2, 5.0)) methods = list(orbital_table["Method"].drop_duplicates()) for x, method in enumerate(methods): subset = orbital_table[orbital_table["Method"] == method] for _, row in subset.iterrows(): label = str(row["Orbital"]) energy = float(row["Energy (eV)"]) frontier = label in {"HOMO", "LUMO"} color = "#2166ac" if float(row["Occupation"]) > 0 else "#b2182b" linewidth = 3.0 if frontier else 1.5 ax.hlines(energy, x - 0.23, x + 0.23, color=color, linewidth=linewidth) ax.text(x + 0.27, energy, label, va="center", fontsize=8) ax.set_xlim(-0.55, len(methods) - 0.25) ax.set_xticks(range(len(methods)), methods) ax.set_ylabel("Orbital energy (eV)") ax.set_title("Frontier orbital energy levels") ax.grid(axis="y", color="#dddddd", linewidth=0.6) fig.tight_layout() fig.savefig(FIGURES / "frontier_energy_levels.png", dpi=200) plt.close(fig) def make_amsreport_orbital(engine_rkf: Path, method_slug: str, orbital: str) -> Path: output_dir = FIGURES / f"amsreport_{method_slug}_{orbital.lower()}" if output_dir.exists(): shutil.rmtree(output_dir) output_dir.mkdir(parents=True) command = [ os.environ["AMSBIN"] + "/amsreport", str(engine_rkf), orbital, "-o", "report.html", "-v", "-scmgeometry 500x400", "-v", "-grid Fine", "-v", "-antialias", "-v", "-bgcolor #ffffff", "-v", "-viewplane {1 2 5}", ] subprocess.run(command, cwd=output_dir, check=True, capture_output=True, text=True) images = sorted(output_dir.rglob("*.jpg")) + sorted(output_dir.rglob("*.png")) if not images: raise RuntimeError(f"amsreport did not create an image for {method_slug} {orbital}") return images[0].relative_to(ROOT) def main() -> None: FIGURES.mkdir(exist_ok=True) TABLES.mkdir(exist_ok=True) workdir = latest_workdir() # Geometry optimization: establishes the one structure used by both single points. optimization = AMSJob.load_external(str(workdir / "ethylene_opt_pbe_dzp")) # PBE single point: supplies the GGA orbital energies and isosurfaces. pbe = AMSJob.load_external(str(workdir / "ethylene_sp_pbe_dzp")) # B3LYP single point: supplies the hybrid orbital energies and isosurfaces. b3lyp = AMSJob.load_external(str(workdir / "ethylene_sp_b3lyp_dzp")) jobs = {"PBE/DZP": pbe, "B3LYP/DZP": b3lyp} summary_rows: list[dict[str, object]] = [] orbital_rows: list[dict[str, object]] = [] for method, job in jobs.items(): homo = float(job.results.get_homo_energies(unit="eV")[0]) lumo = float(job.results.get_lumo_energies(unit="eV")[0]) gap = float(job.results.get_smallest_homo_lumo_gap(unit="eV")) summary_rows.append( {"Method": method, "HOMO (eV)": homo, "LUMO (eV)": lumo, "Gap (eV)": gap} ) orbital_rows.extend(frontier_rows(job, method)) summary = pd.DataFrame(summary_rows) orbitals = pd.DataFrame(orbital_rows) summary.to_csv(TABLES / "frontier_summary.csv", index=False) orbitals.to_csv(TABLES / "frontier_orbitals.csv", index=False) make_energy_level_plot(orbitals) optimized_system = optimization.results.get_main_system() view( optimized_system, guess_bonds=len(optimized_system.bonds) == 0, direction="along_pca3", width=600, height=450, picture_path=str(FIGURES / "optimized_ethylene.png"), ) orbital_images: dict[tuple[str, str], Path] = {} for method_slug, job in (("pbe", pbe), ("b3lyp", b3lyp)): engine_rkf = Path(job.results.rkfpath(file="engine")).resolve() for orbital in ("HOMO", "LUMO"): orbital_images[(method_slug, orbital)] = make_amsreport_orbital( engine_rkf, method_slug, orbital ) carbon_indices = [i for i, atom in enumerate(optimized_system.atoms) if atom.symbol == "C"] cc_distance = optimized_system.get_distance(carbon_indices[0], carbon_indices[1]) gap_change = float(summary.loc[summary["Method"] == "B3LYP/DZP", "Gap (eV)"].iloc[0]) - float( summary.loc[summary["Method"] == "PBE/DZP", "Gap (eV)"].iloc[0] ) lines = [ "# Ethylene frontier orbitals with PBE and B3LYP", "", "Ethylene was built from the SMILES string `C=C` and optimized once with ADF at PBE/DZP. " "Two single-point calculations, PBE/DZP and B3LYP/DZP, then used that same geometry. " "Both calculations used an all-electron DZP basis and Normal numerical quality. Orbital energies are Kohn-Sham eigenvalues in eV.", "", "## Optimized geometry", "", f"The optimized C=C distance is {cc_distance:.4f} Å.", "", "![Optimized ethylene geometry](figures/optimized_ethylene.png)", "", "## Frontier orbital energies", "", summary.to_markdown(index=False, floatfmt=".4f"), "", "The nearest three occupied and three virtual orbitals are:", "", orbitals.to_markdown(index=False, floatfmt=("", "", ".4f", ".1f")), "", "![Frontier orbital energy-level diagram](figures/frontier_energy_levels.png)", "", "## HOMO and LUMO isosurfaces", "", "The four orbital pictures below were generated by `amsreport` from the corresponding ADF engine result files. Contrasting colors mark opposite orbital phases.", "", "| Method | HOMO | LUMO |", "|---|---|---|", f"| PBE/DZP | ![PBE HOMO]({orbital_images[('pbe', 'HOMO')].as_posix()}) | ![PBE LUMO]({orbital_images[('pbe', 'LUMO')].as_posix()}) |", f"| B3LYP/DZP | ![B3LYP HOMO]({orbital_images[('b3lyp', 'HOMO')].as_posix()}) | ![B3LYP LUMO]({orbital_images[('b3lyp', 'LUMO')].as_posix()}) |", "", "## Conclusion", "", f"At this fixed PBE geometry, replacing PBE with B3LYP changes the Kohn-Sham HOMO-LUMO gap by {gap_change:+.4f} eV. " "The HOMO and LUMO isosurfaces allow the comparison to distinguish an energy shift from a qualitative change in orbital shape. " "These Kohn-Sham gaps are method-dependent orbital-energy differences, not optical excitation energies.", "", "## Calculation inputs", "", ] for title, job in ( ("PBE/DZP geometry optimization", optimization), ("PBE/DZP single point", pbe), ("B3LYP/DZP single point", b3lyp), ): lines.extend([f"### {title}", "", "```ams", job.get_input().rstrip(), "```", ""]) (ROOT / "report.md").write_text("\n".join(lines), encoding="utf-8") if __name__ == "__main__": main()