#!/usr/bin/env amspython from __future__ import annotations from pathlib import Path import matplotlib.pyplot as plt import numpy as np import pandas as pd from scm.base import Units from scm.plams import AMSJob, view ROOT = Path(__file__).resolve().parent SMILES = "C=CC(CCCC)C1C=CC(C)=CC1=O" def latest_workdir(stem: str) -> Path: candidates = [path for path in ROOT.glob(f"{stem}*") if path.is_dir()] if not candidates: raise FileNotFoundError(f"No work directory matching {stem!r}") def generation(path: Path) -> int: suffix = path.name[len(stem) :] return 1 if not suffix else int(suffix.lstrip(".")) return max(candidates, key=generation) def gibbs_hartree(job: AMSJob) -> float: return float(job.results.readrkf("Thermodynamics", "Gibbs free Energy", file="engine")) def render_system(system: object, path: Path, width: int = 480, height: int = 360) -> None: view( system, guess_bonds=len(system.bonds) == 0, direction="along_pca3", width=width, height=height, picture_path=str(path), ) def main() -> None: workdir = latest_workdir("01-run_workdir") figures = ROOT / "figures" tables = ROOT / "tables" figures.mkdir(exist_ok=True) tables.mkdir(exist_ok=True) # Each job is loaded directly from its AMS results directory. The accompanying # string documents why that calculation supports the report. jobs_with_justification = { "Reactant": ( AMSJob.load_external(str(workdir / "reactant_opt_freq")), "GFN1-xTB reactant geometry optimization and harmonic normal modes.", ), "TS": ( AMSJob.load_external(str(workdir / "ts_opt_freq")), "GFN1-xTB transition-state search from the highest scan-energy geometry, with calculated initial Hessian and final normal modes.", ), "Product": ( AMSJob.load_external(str(workdir / "product_opt_freq")), "GFN1-xTB geometry optimization and harmonic normal modes from the 1.5 Å scan endpoint.", ), "Scan": ( AMSJob.load_external(str(workdir / "oc_distance_scan")), "Relaxed O···C distance scan used to locate the TS-search starting geometry.", ), } reactant_job = jobs_with_justification["Reactant"][0] ts_job = jobs_with_justification["TS"][0] product_job = jobs_with_justification["Product"][0] scan_job = jobs_with_justification["Scan"][0] hartree_to_kcal = Units.conversion_factor("hartree", "kcal/mol") bohr_to_angstrom = Units.conversion_factor("bohr", "angstrom") rows: list[dict[str, object]] = [] for label, job in (("Reactant", reactant_job), ("TS", ts_job), ("Product", product_job)): frequencies = np.asarray(job.results.get_frequencies(unit="cm^-1"), dtype=float) rows.append( { "Structure": label, "Electronic energy (Hartree)": float(job.results.get_energy(unit="hartree")), "Gibbs free energy (Hartree)": gibbs_hartree(job), "PES point character": str( job.results.readrkf("AMSResults", "PESPointCharacter", file="engine") ), "Negative frequencies": int(np.sum(frequencies < 0.0)), "Lowest frequency (cm^-1)": float(np.min(frequencies)), } ) energy_df = pd.DataFrame(rows) energy_df["Electronic energy (kcal/mol)"] = energy_df["Electronic energy (Hartree)"] * hartree_to_kcal energy_df["Relative electronic energy (kcal/mol)"] = ( energy_df["Electronic energy (Hartree)"] - energy_df.loc[0, "Electronic energy (Hartree)"] ) * hartree_to_kcal energy_df["Gibbs free energy (kcal/mol)"] = energy_df["Gibbs free energy (Hartree)"] * hartree_to_kcal energy_df["Relative Gibbs free energy (kcal/mol)"] = ( energy_df["Gibbs free energy (Hartree)"] - energy_df.loc[0, "Gibbs free energy (Hartree)"] ) * hartree_to_kcal column_order = [ "Structure", "Electronic energy (Hartree)", "Electronic energy (kcal/mol)", "Relative electronic energy (kcal/mol)", "Gibbs free energy (Hartree)", "Gibbs free energy (kcal/mol)", "Relative Gibbs free energy (kcal/mol)", "PES point character", "Negative frequencies", "Lowest frequency (cm^-1)", ] energy_df = energy_df[column_order] energy_df.to_csv(tables / "stationary_points.csv", index=False) scan = scan_job.results.get_pesscan_results() distances = np.asarray(scan["RaveledPESCoords"][0], dtype=float) * bohr_to_angstrom scan_energies = np.asarray(scan["PES"], dtype=float) relative_scan_energies = (scan_energies - scan_energies[0]) * hartree_to_kcal scan_df = pd.DataFrame( { "Point": np.arange(1, len(distances) + 1), "O-C distance (angstrom)": distances, "Electronic energy (Hartree)": scan_energies, "Relative electronic energy (kcal/mol)": relative_scan_energies, "Converged": np.asarray(scan["Converged"], dtype=bool), } ) scan_df.to_csv(tables / "distance_scan.csv", index=False) highest_index = int(np.argmax(scan_energies)) fig, ax = plt.subplots(figsize=(6.4, 4.2)) ax.plot(distances, relative_scan_energies, "o-", color="#2457a6", linewidth=1.5, markersize=4) ax.scatter( [distances[highest_index]], [relative_scan_energies[highest_index]], color="#b2292e", marker="*", s=130, zorder=3, label=f"TS-search start (point {highest_index + 1})", ) ax.set_xlabel("Carbonyl O–terminal vinyl C distance (Å)") ax.set_ylabel("Electronic energy relative to first scan point (kcal/mol)") ax.invert_xaxis() ax.grid(alpha=0.25) ax.legend(frameon=False) fig.tight_layout() fig.savefig(figures / "energy_vs_distance.png", dpi=220) plt.close(fig) stationary_images: list[tuple[str, Path]] = [] for label, job in (("Reactant", reactant_job), ("TS", ts_job), ("Product", product_job)): image_path = figures / f"{label.lower()}_optimized.png" system = job.results.get_main_system() render_system(system, image_path) stationary_images.append((label, image_path)) representative_indices = np.linspace(0, len(distances) - 1, 5).round().astype(int).tolist() representative_paths: list[Path] = [] for sequence, point_index in enumerate(representative_indices, start=1): molecule = scan["Molecules"][point_index] image_path = figures / f"scan_representative_{sequence}.png" render_system(molecule, image_path, width=360, height=300) representative_paths.append(image_path) fig, axes = plt.subplots(1, 5, figsize=(15, 3.2)) for ax, image_path, point_index in zip(axes, representative_paths, representative_indices): ax.imshow(plt.imread(image_path)) ax.set_title(f"{distances[point_index]:.1f} Å\npoint {point_index + 1}", fontsize=10) ax.axis("off") fig.tight_layout(w_pad=0.1) fig.savefig(figures / "scan_structures_one_row.png", dpi=180, bbox_inches="tight") plt.close(fig) barrier = float(energy_df.loc[energy_df["Structure"] == "TS", "Relative Gibbs free energy (kcal/mol)"].iloc[0]) reaction_free_energy = float( energy_df.loc[energy_df["Structure"] == "Product", "Relative Gibbs free energy (kcal/mol)"].iloc[0] ) reverse_barrier = barrier - reaction_free_energy reverse_reaction_free_energy = -reaction_free_energy ts_frequencies = np.asarray(ts_job.results.get_frequencies(unit="cm^-1"), dtype=float) ts_negative = ts_frequencies[ts_frequencies < 0.0] reactant_frequencies = np.asarray(reactant_job.results.get_frequencies(unit="cm^-1"), dtype=float) product_frequencies = np.asarray(product_job.results.get_frequencies(unit="cm^-1"), dtype=float) display_energy_df = energy_df.copy() numeric_columns = display_energy_df.select_dtypes(include=[np.number]).columns display_energy_df[numeric_columns] = display_energy_df[numeric_columns].round(4) display_scan_df = scan_df.copy() display_scan_df["O-C distance (angstrom)"] = display_scan_df["O-C distance (angstrom)"].round(4) display_scan_df["Electronic energy (Hartree)"] = display_scan_df["Electronic energy (Hartree)"].round(8) display_scan_df["Relative electronic energy (kcal/mol)"] = display_scan_df[ "Relative electronic energy (kcal/mol)" ].round(4) inputs = [] for label in ("Reactant", "Scan", "TS", "Product"): job, justification = jobs_with_justification[label] inputs.append( f"### {label}\n\n{justification}\n\n```ams\n{job.get_input().rstrip()}\n```" ) report = f"""# GFN1-xTB Claisen rearrangement profile ## Summary The neutral, closed-shell molecule `{SMILES}` was treated in the gas phase with GFN1-xTB. The carbonyl oxygen and terminal vinyl carbon were identified from the `ChemicalSystem.bonds` graph as atoms 15 and 1, respectively. All thermochemical values are ideal-gas rigid-rotor/harmonic-oscillator results at the AMS default temperature of 298.15 K. For the explicitly scanned direction—supplied dienone → short O–C aryl-ether endpoint—the reaction free-energy barrier is **{barrier:.2f} kcal/mol**, and the reaction free energy is **{reaction_free_energy:.2f} kcal/mol**. The requested O–C contraction is the retro-Claisen direction. For the opposite, conventional aryl ether → dienone Claisen direction, the same stationary points give ΔG‡ = **{reverse_barrier:.2f} kcal/mol** and ΔG = **{reverse_reaction_free_energy:.2f} kcal/mol**. The transition state has exactly one imaginary frequency, **{ts_negative[0]:.2f} cm⁻¹**, and AMS classifies it as a transition state. ## Stationary points | Reactant | Transition state | Product | |:--:|:--:|:--:| | ![Optimized reactant](figures/reactant_optimized.png) | ![Optimized transition state](figures/ts_optimized.png) | ![Optimized product](figures/product_optimized.png) | {display_energy_df.to_markdown(index=False)} Electronic and Gibbs energies use identical GFN1-xTB settings. Relative values are referenced to the optimized reactant. The reactant has no negative frequencies and is classified as a local minimum. The product is classified as a local minimum by AMS; its lowest mode is {np.min(product_frequencies):.2f} cm⁻¹. The small negative value is below the AMS PES-point-character tolerance and represents a very soft conformational mode rather than a chemically significant unstable mode. ## Relaxed O–C distance scan The optimized reactant O···C distance is 4.758 Å. To preserve an exact 0.2 Å grid ending at 1.5 Å, the relaxed scan begins at 4.9 Å, requiring a 0.142 Å adjustment at its first constrained point. The only scan warning records this difference between the initial geometry and first requested constraint. All 18 constrained optimizations converged. Point {highest_index + 1}, at {distances[highest_index]:.2f} Å, has the highest scan energy and was used directly as the TS-search starting structure. ![Electronic energy versus O-C distance](figures/energy_vs_distance.png) Five representative relaxed-scan geometries, displayed in scan order in one row: ![Five representative three-dimensional scan structures](figures/scan_structures_one_row.png) {display_scan_df.to_markdown(index=False)} ## Frequency and TS verification - Reactant: {len(reactant_frequencies)} vibrational modes; {int(np.sum(reactant_frequencies < 0.0))} negative frequencies; AMS character `local minimum`. - Transition state: {len(ts_frequencies)} vibrational modes; {len(ts_negative)} negative frequency ({ts_negative[0]:.2f} cm⁻¹); AMS character `transition state`. - Product: {len(product_frequencies)} vibrational modes; {int(np.sum(product_frequencies < 0.0))} formally negative frequency ({np.min(product_frequencies):.2f} cm⁻¹); AMS character `local minimum`. The TS calculation explicitly used `GeometryOptimization InitialHessian Type Calculate`, requested final normal modes, and followed the O–C distance reaction coordinate. ## Conclusion At the GFN1-xTB/RRHO level, the requested O–C contraction has ΔG‡ = **{barrier:.2f} kcal/mol** and ΔG = **{reaction_free_energy:.2f} kcal/mol** at 298.15 K; the opposite conventional Claisen direction has ΔG‡ = **{reverse_barrier:.2f} kcal/mol** and ΔG = **{reverse_reaction_free_energy:.2f} kcal/mol**. The stationary-point classification and one-imaginary-frequency requirement are satisfied for the transition state. Because the molecule is conformationally flexible and only one UFF-generated starting conformer was refined, these values are a single-conformer estimate. ## Calculation inputs and provenance {chr(10).join(inputs)} """ (ROOT / "report.md").write_text(report) print(f"Wrote {ROOT / 'report.md'}") print(f"Barrier: {barrier:.4f} kcal/mol") print(f"Reaction free energy: {reaction_free_energy:.4f} kcal/mol") print(f"TS negative frequencies: {ts_negative.tolist()}") if __name__ == "__main__": main()