#!/usr/bin/env amspython from __future__ import annotations from pathlib import Path from typing import Dict, Optional import numpy as np from scm.plams import AMSJob, AMSMSDJob, AMSNVTJob, Settings, finish, init, log from common import ( SAMPLING_FREQUENCY, TARGET_TEMPERATURE_K, TIMESTEP_FS, completed_jobs, fragmentation_failure_paths, restart_settings, sorted_cumulative_jobs, ) INITIAL_PS = 100 INCREMENT_PS = 50 MAXIMUM_PS = 1000 MSD_WINDOW_FS = 20_000.0 MSD_FIT_START_FS = 5_000.0 def start_production(equilibration_job: AMSJob) -> AMSJob: overrides = Settings() job = AMSNVTJob.restart_from( equilibration_job, name="prod_100ps", settings=overrides, nsteps=int(INITIAL_PS * 1000 / TIMESTEP_FS), timestep=TIMESTEP_FS, samplingfreq=SAMPLING_FREQUENCY, temperature=TARGET_TEMPERATURE_K, thermostat="NHC", tau=100.0, calcpressure=True, writevelocities=True, writebonds=True, writemolecules=True, ) md = job.settings.input.ams.MolecularDynamics if "Restart" in md: del md.Restart if "CopyRestartTrajectory" in md: del md.CopyRestartTrajectory if "LoadSystem" in job.settings.input.ams: del job.settings.input.ams.LoadSystem log("Starting a separate 100 ps production NVT trajectory with an NHC thermostat. " "Only the final equilibration coordinates and velocities are loaded; equilibration frames are not copied.") job.run() if not job.ok(): raise RuntimeError(f"Initial production trajectory failed: {job.results.get_errormsg()}") log(f"Completed the cumulative 100 ps production trajectory at {job.path}.") return job def extend(previous_job: AMSJob, target_ps: int) -> AMSJob: settings = restart_settings(previous_job, int(target_ps * 1000 / TIMESTEP_FS)) job = AMSJob(name=f"prod_{target_ps:04d}ps", molecule=None, settings=settings) log(f"Extending production exactly from {previous_job.name} to {target_ps} ps with Restart and " "CopyRestartTrajectory Yes.") job.run() if not job.ok(): raise RuntimeError(f"Production extension to {target_ps} ps failed: {job.results.get_errormsg()}") log(f"Completed the cumulative {target_ps} ps production trajectory at {job.path}.") return job def analysis_settings() -> Settings: settings = Settings() msd = settings.input.MeanSquareDisplacement msd.Atoms.Region = "glycine" msd.UseMolecularCentersOfMass.Enabled = "Yes" msd.UseMolecularCentersOfMass.CheckForBondChanges = "Yes" msd.ComputeStandardDeviation = "Yes" return settings def reusable_msd_job(name: str) -> Optional[AMSMSDJob]: candidates = [] from common import ROOT, natural_key for workdir in ROOT.glob("03-production_workdir*"): path = workdir / name if path.is_dir() and (path / f"{name}.err").is_file(): candidates.append(path) for path in sorted(candidates, key=lambda item: natural_key(str(item))): try: job = AMSMSDJob.load_external(str(path)) if job.ok(): return job except Exception as exc: log(f"Could not reuse MSD analysis at {path}: {exc}") return None def run_msd(production_job: AMSJob, duration_ps: int) -> Optional[float]: name = f"msd_{duration_ps:04d}ps" job = reusable_msd_job(name) if job is not None: value = float(job.results.get_diffusion_coefficient( start_time_fit_fs=MSD_FIT_START_FS, end_time_fit_fs=MSD_WINDOW_FS )) log(f"Reusing {name}: D = {value:.6e} m^2/s from the 5-20 ps fit interval.") return value log(f"Running glycine molecular-center-of-mass MSD analysis for the cumulative {duration_ps} ps trajectory. " "The correlation window is 20 ps and the linear fit interval is 5-20 ps.") job = AMSMSDJob( previous_job=production_job, name=name, max_correlation_time_fs=MSD_WINDOW_FS, start_time_fit_fs=MSD_FIT_START_FS, settings=analysis_settings(), ) try: job.run() if not job.ok(): text = " ".join(job.results.grep_output("bond", options="-i") or []) if "bond" in text.lower(): log("Molecular-center-of-mass MSD failed because the glycine bonding changed. " "Production will stop and the fragmentation status will be reported.") return None raise RuntimeError(f"MSD analysis failed: {job.results.get_errormsg()}") value = float(job.results.get_diffusion_coefficient( start_time_fit_fs=MSD_FIT_START_FS, end_time_fit_fs=MSD_WINDOW_FS )) except Exception as exc: text = str(exc).lower() job_text = "" if job.path: from common import read_job_text job_text = read_job_text(Path(job.path)).lower() combined = text + "\n" + job_text if "bond" in combined and any(word in combined for word in ("change", "changed", "molecule", "fragment")): log(f"Molecular-center-of-mass MSD detected glycine fragmentation or atom exchange: {exc}") return None raise log(f"Cumulative {duration_ps} ps diffusion estimate: D = {value:.6e} m^2/s.") return value def converged(estimates: Dict[int, float]) -> bool: if len(estimates) < 3: log(f"Production convergence needs three cumulative estimates; only {len(estimates)} are available.") return False durations = sorted(estimates)[-3:] values = np.asarray([estimates[duration] for duration in durations], dtype=float) mean = float(values.mean()) relative = np.abs(values - mean) / abs(mean) if mean != 0.0 else np.full(3, np.inf) for duration, value, deviation in zip(durations, values, relative): log(f"Production convergence member at {duration} ps: D = {value:.6e} m^2/s; " f"deviation from three-estimate mean = {100 * deviation:.3f}%; required <= 15.000%.") passed = bool(np.all(relative <= 0.15)) log(f"Production convergence check using {durations}: mean D = {mean:.6e} m^2/s; " f"{'PASS' if passed else 'FAIL'}.") return passed def main() -> None: prior_fragmentation = fragmentation_failure_paths() if prior_fragmentation: log(f"A prior center-of-mass analysis recorded bond changes at {prior_fragmentation[-1]}. " "No further production MD will be run.") return equilibration = sorted_cumulative_jobs(completed_jobs("02-equilibration_workdir", "eq_*ps"), "eq_") if not equilibration: raise RuntimeError("No completed equilibration job. Run 02-equilibration.py first.") final_equilibration = equilibration[-1][1] log(f"Using the final cumulative equilibration trajectory {final_equilibration.name} at " f"{final_equilibration.path} to initialize production.") cumulative = sorted_cumulative_jobs(completed_jobs("03-production_workdir", "prod_*ps"), "prod_") if cumulative: log(f"Found {len(cumulative)} reusable completed production trajectories; latest is " f"{cumulative[-1][0]} ps at {cumulative[-1][1].path}.") else: first = start_production(final_equilibration) cumulative = [(INITIAL_PS, first)] estimates: Dict[int, float] = {} for duration_ps, job in cumulative: value = run_msd(job, duration_ps) if value is None: return estimates[duration_ps] = value if converged(estimates): log(f"Production already converged at {max(estimates)} ps; no extension is needed.") return while cumulative[-1][0] < MAXIMUM_PS: target_ps = min(cumulative[-1][0] + INCREMENT_PS, MAXIMUM_PS) job = extend(cumulative[-1][1], target_ps) cumulative.append((target_ps, job)) value = run_msd(job, target_ps) if value is None: return estimates[target_ps] = value if converged(estimates): log(f"Production diffusion estimate converged at {target_ps} ps.") return log(f"Production reached the 1 ns cap without convergence. The latest estimate, " f"D = {estimates[MAXIMUM_PS]:.6e} m^2/s, will be reported as unconverged.") if __name__ == "__main__": init(folder="03-production_workdir") try: main() finally: finish()