#!/usr/bin/env amspython from __future__ import annotations from pathlib import Path from typing import Dict, List, Optional, Tuple import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import pandas as pd from scm.base import Units from scm.plams import AMSJob, AMSMSDJob, log, view from common import ( ROOT, TARGET_TEMPERATURE_K, completed_jobs, equilibration_metrics, fragmentation_failure_paths, natural_key, sorted_cumulative_jobs, ) TABLES = ROOT / "tables" FIGURES = ROOT / "figures" MSD_FIT_START_FS = 5_000.0 MSD_WINDOW_FS = 20_000.0 def successful_msd_jobs() -> Dict[int, AMSMSDJob]: found: Dict[int, AMSMSDJob] = {} candidates: List[Path] = [] for workdir in ROOT.glob("03-production_workdir*"): candidates.extend(path for path in workdir.glob("msd_*ps") if path.is_dir()) for path in sorted(candidates, key=lambda item: natural_key(str(item))): try: job = AMSMSDJob.load_external(str(path)) if job.ok(): duration = int(path.name.split("_")[1][:-2]) found[duration] = job except Exception as exc: log(f"Skipping unreadable MSD analysis at {path}: {exc}") return found def equilibration_table(equilibration: List[Tuple[int, AMSJob]]) -> pd.DataFrame: rows: List[dict] = [] streak = 0 for duration, job in equilibration: if duration < 10: continue metrics = equilibration_metrics(job) streak = streak + 1 if metrics["passed"] else 0 rows.append({ "duration_ps": duration, "window_start_ps": float(metrics["window_start_fs"]) / 1000, "window_end_ps": float(metrics["window_end_fs"]) / 1000, "fitted_energy_change_hartree": metrics["fitted_change_hartree"], "energy_std_hartree": metrics["trajectory_std_hartree"], "abs_change_over_std": metrics["energy_change_in_std"], "mean_temperature_k": metrics["mean_temperature_k"], "temperature_error_percent": 100 * float(metrics["temperature_relative_error"]), "energy_pass": metrics["energy_ok"], "temperature_pass": metrics["temperature_ok"], "combined_pass": metrics["passed"], "consecutive_passes": streak, "accepted": streak >= 2, }) return pd.DataFrame(rows) def production_table(msd_jobs: Dict[int, AMSMSDJob]) -> pd.DataFrame: rows: List[dict] = [] estimates: Dict[int, float] = {} for duration in sorted(msd_jobs): value = float(msd_jobs[duration].results.get_diffusion_coefficient( start_time_fit_fs=MSD_FIT_START_FS, end_time_fit_fs=MSD_WINDOW_FS )) estimates[duration] = value last = sorted(estimates)[-3:] mean: Optional[float] = None maximum_deviation: Optional[float] = None converged = False if len(last) == 3: values = np.asarray([estimates[item] for item in last]) mean = float(values.mean()) deviations = np.abs(values - mean) / abs(mean) if mean != 0.0 else np.full(3, np.inf) maximum_deviation = float(deviations.max()) converged = bool(np.all(deviations <= 0.15)) rows.append({ "duration_ps": duration, "diffusion_m2_s": value, "last_three_mean_m2_s": mean, "last_three_max_deviation_percent": None if maximum_deviation is None else 100 * maximum_deviation, "converged": converged, }) return pd.DataFrame(rows) def history_frame(job: AMSJob, stage: str) -> pd.DataFrame: time = np.asarray(job.results.get_history_property("Time", "MDHistory"), dtype=float) energy = np.asarray(job.results.get_history_property("PotentialEnergy", "MDHistory"), dtype=float) temperature = np.asarray(job.results.get_history_property("Temperature", "MDHistory"), dtype=float) pressure_au = np.asarray(job.results.get_history_property("Pressure", "MDHistory"), dtype=float) energy_factor = Units.conversion_factor("hartree", "kJ/mol") pressure_factor = Units.conversion_factor("au", "Pa") return pd.DataFrame({ "stage": stage, "time_fs": time, "potential_energy_hartree": energy, "potential_energy_kj_mol": energy * energy_factor, "temperature_k": temperature, "pressure_pa": pressure_au * pressure_factor, }) def make_trace_figure(eq_trace: pd.DataFrame, prod_trace: pd.DataFrame) -> None: fig, axes = plt.subplots(2, 1, figsize=(8, 7), sharex=False) for frame, label in ((eq_trace, "Equilibration"), (prod_trace, "Production")): axes[0].plot(frame["time_fs"] / 1000, frame["potential_energy_kj_mol"], label=label) axes[1].plot(frame["time_fs"] / 1000, frame["temperature_k"], label=label) axes[0].set_ylabel("Potential energy (kJ/mol)") axes[1].set_ylabel("Temperature (K)") axes[1].set_xlabel("Trajectory time (ps)") axes[1].axhline(TARGET_TEMPERATURE_K, color="black", linewidth=0.8, linestyle=":") for axis in axes: axis.legend() axis.grid(alpha=0.2) fig.tight_layout() fig.savefig(FIGURES / "potential-energy-temperature.png", dpi=180) plt.close(fig) def make_msd_figure(job: AMSMSDJob) -> pd.DataFrame: time, msd = job.results.get_msd() fit, fit_time, fit_msd = job.results.get_linear_fit( start_time_fit_fs=MSD_FIT_START_FS, end_time_fit_fs=MSD_WINDOW_FS ) frame = pd.DataFrame({"time_fs": time, "msd_angstrom2": msd}) frame.to_csv(TABLES / "glycine-com-msd.csv", index=False) fig, ax = plt.subplots(figsize=(7, 4.5)) ax.plot(time / 1000, msd, label="Glycine COM MSD") ax.plot(fit_time / 1000, fit_msd, label="Linear fit, 5-20 ps") ax.axvspan(5, 20, color="tab:orange", alpha=0.12) ax.set_xlabel("Correlation time (ps)") ax.set_ylabel("MSD (angstrom^2)") ax.grid(alpha=0.2) ax.legend() fig.tight_layout() fig.savefig(FIGURES / "glycine-com-msd-fit.png", dpi=180) plt.close(fig) return frame def input_section(title: str, job: AMSJob) -> str: return f"### {title}\n\nProvenance: `{job.path}`\n\n```ams\n{job.get_input().rstrip()}\n```\n" def main() -> None: TABLES.mkdir(exist_ok=True) FIGURES.mkdir(exist_ok=True) setup_map = completed_jobs("01-initial-system_workdir", "setup_opt") equilibration = sorted_cumulative_jobs(completed_jobs("02-equilibration_workdir", "eq_*ps"), "eq_") production = sorted_cumulative_jobs(completed_jobs("03-production_workdir", "prod_*ps"), "prod_") msd_jobs = successful_msd_jobs() fragmentation = fragmentation_failure_paths() if not setup_map or not equilibration or not production: raise RuntimeError("Setup, equilibration, and production jobs must exist before rendering the report.") setup = setup_map["setup_opt"] final_eq_duration, final_eq = equilibration[-1] final_prod_duration, final_prod = production[-1] eq_table = equilibration_table(equilibration) prod_table = production_table(msd_jobs) eq_table.to_csv(TABLES / "equilibration-convergence.csv", index=False) prod_table.to_csv(TABLES / "production-convergence.csv", index=False) eq_trace = history_frame(final_eq, "equilibration") prod_trace = history_frame(final_prod, "production") eq_trace.to_csv(TABLES / "equilibration-trace.csv", index=False) prod_trace.to_csv(TABLES / "production-trace.csv", index=False) make_trace_figure(eq_trace, prod_trace) packed_system = setup.results.get_input_system() view( packed_system, guess_bonds=len(packed_system.bonds) == 0, direction="tilt_x", width=700, height=700, picture_path=str(FIGURES / "packed-system.png"), ) final_msd: Optional[AMSMSDJob] = None diffusion: Optional[float] = None fit_r_squared: Optional[float] = None if msd_jobs: latest_msd_duration = max(msd_jobs) final_msd = msd_jobs[latest_msd_duration] diffusion = float(final_msd.results.get_diffusion_coefficient( start_time_fit_fs=MSD_FIT_START_FS, end_time_fit_fs=MSD_WINDOW_FS )) fit_result, _, _ = final_msd.results.get_linear_fit( start_time_fit_fs=MSD_FIT_START_FS, end_time_fit_fs=MSD_WINDOW_FS ) fit_r_squared = float(fit_result.rvalue ** 2) make_msd_figure(final_msd) eq_converged = bool(len(eq_table) and eq_table.iloc[-1]["accepted"]) prod_converged = bool(len(prod_table) and prod_table.iloc[-1]["converged"]) fragmentation_status = "detected" if fragmentation else "not detected" prod_temp_mean = float(prod_trace["temperature_k"].mean()) prod_temp_std = float(prod_trace["temperature_k"].std(ddof=1)) prod_pressure_mean = float(prod_trace["pressure_pa"].mean()) prod_pressure_std = float(prod_trace["pressure_pa"].std(ddof=1)) density = setup.results.get_input_molecule().get_density() / 1000 if diffusion is None: result_text = "No valid molecular-center-of-mass diffusion estimate was produced." else: result_text = (f"The final glycine center-of-mass diffusion coefficient is **{diffusion:.6e} m^2/s** " f"with a fitted MSD R^2 of {fit_r_squared:.6f}.") if fragmentation: convergence_text = "The molecular-center-of-mass analysis detected a glycine bond change, so production stopped and the result is not converged." elif prod_converged: convergence_text = "The last three cumulative estimates are each within 15% of their mean, so the production estimate is converged by the specified test." elif final_prod_duration >= 1000: convergence_text = "The trajectory reached the 1 ns cap without meeting the three-estimate test. The latest estimate is reported as unconverged." else: convergence_text = "The available trajectory has not yet met the three-estimate convergence test." markdown = f"""# Glycine diffusion in water at 380 K ## Result {result_text} The MSD correlation window is 20 ps and the fitted interval is 5-20 ps. The cumulative production trajectory is {final_prod_duration} ps. {convergence_text} The mean of the final three cumulative estimates is {float(prod_table.iloc[-1]['last_three_mean_m2_s']):.6e} m^2/s. This mean is used only for the convergence test; the reported coefficient comes from the full {final_prod_duration} ps trajectory. The production trajectory has a mean temperature of {prod_temp_mean:.3f} +/- {prod_temp_std:.3f} K and a mean pressure of {prod_pressure_mean:.6e} +/- {prod_pressure_std:.6e} Pa. The fragmentation or atom-exchange status is **{fragmentation_status}**. ## System and protocol The periodic cubic starting system contains one neutral NH2-CH2-COOH glycine molecule and 40 water molecules. It has {len(packed_system)} atoms and an initial packed density of {density:.6f} g/cm^3. ReaxFF used `Glycine.ff` with `nproc=1`. The geometry optimization used Basic convergence. Equilibration used Berendsen NVT at 380 K with a 100 fs coupling time. Production used a separate NHC NVT trajectory with a 100 fs coupling time. Both MD stages used a 0.5 fs timestep and stored a frame every 100 fs. The equilibration endpoint is {final_eq_duration} ps and its convergence status is **{'accepted' if eq_converged else 'not accepted by the two-check test'}**. ![Packed periodic system](figures/packed-system.png) ## Equilibration convergence {eq_table.to_markdown(index=False)} The energy test fits the five block-mean potential energies over the trailing half of each accumulated trajectory. It passes when the fitted change is no more than 0.5 times the trajectory standard deviation. The temperature test requires the trailing mean to lie within 2% of 380 K. Two consecutive combined passes are required. ## Production convergence {prod_table.to_markdown(index=False) if len(prod_table) else 'No successful center-of-mass MSD analyses were available.'} ![Potential-energy and temperature traces](figures/potential-energy-temperature.png) {'![Glycine center-of-mass MSD and fitted region](figures/glycine-com-msd-fit.png)' if final_msd is not None else ''} CSV data are in [`tables/`](tables/). The convergence histories, thermodynamic traces, and MSD values shown here come directly from the individual AMS and Analysis job directories. ## Limitations The box contains only 40 water molecules. Periodic hydrodynamic finite-size effects can therefore bias the diffusion coefficient. No finite-size correction or independent replicate trajectories were applied. The fixed-density NVT state also prevents density relaxation and makes the reported mean pressure an observation rather than a controlled state variable. ## Complete calculation inputs """ markdown += input_section("Initial ReaxFF geometry optimization", setup) for duration, job in equilibration: markdown += input_section(f"Equilibration, cumulative {duration} ps", job) for duration, job in production: markdown += input_section(f"Production, cumulative {duration} ps", job) if final_msd is not None: input_files = sorted(Path(final_msd.path).glob("*.in")) if input_files: markdown += f"### Final MSD analysis\n\nProvenance: `{final_msd.path}`\n\n```ams\n{input_files[0].read_text().rstrip()}\n```\n" markdown += "\n## Conclusion\n\n" + result_text + " " + convergence_text + "\n" (ROOT / "report.md").write_text(markdown) log(f"Wrote {(ROOT / 'report.md')} with tables and figures.") if __name__ == "__main__": main()