from __future__ import annotations import re from pathlib import Path from typing import Dict, List, Optional, Tuple import numpy as np from scm.plams import AMSJob, Settings, log ROOT = Path(__file__).resolve().parent FORCE_FIELD = "Glycine.ff" TIMESTEP_FS = 0.5 SAMPLING_INTERVAL_FS = 100.0 SAMPLING_FREQUENCY = int(SAMPLING_INTERVAL_FS / TIMESTEP_FS) TARGET_TEMPERATURE_K = 380.0 def natural_key(value: str) -> List[object]: return [int(part) if part.isdigit() else part.lower() for part in re.split(r"(\d+)", value)] def job_paths(workdir_prefix: str, name_pattern: str) -> List[Path]: paths: List[Path] = [] for workdir in ROOT.glob(f"{workdir_prefix}*"): if not workdir.is_dir(): continue for path in workdir.glob(name_pattern): if path.is_dir() and (path / "ams.rkf").is_file(): paths.append(path) return sorted(paths, key=lambda path: natural_key(str(path))) def load_completed(path: Path) -> Optional[AMSJob]: try: job = AMSJob.load_external(str(path)) if job.ok(): return job except Exception as exc: log(f"Could not reuse {path}: {exc}") return None def completed_jobs(workdir_prefix: str, name_pattern: str) -> Dict[str, AMSJob]: found: Dict[str, AMSJob] = {} for path in job_paths(workdir_prefix, name_pattern): job = load_completed(path) if job is not None: found[path.name] = job return found def reaxff_settings() -> Settings: settings = Settings() settings.input.ReaxFF.ForceField = FORCE_FIELD settings.runscript.nproc = 1 return settings def duration_from_name(name: str) -> int: match = re.search(r"_(\d+)ps$", name) if match is None: raise ValueError(f"Cannot read trajectory duration from job name {name!r}") return int(match.group(1)) def sorted_cumulative_jobs(jobs: Dict[str, AMSJob], prefix: str) -> List[Tuple[int, AMSJob]]: items = [(duration_from_name(name), job) for name, job in jobs.items() if name.startswith(prefix)] return sorted(items, key=lambda item: item[0]) def equilibration_metrics(job: AMSJob) -> Dict[str, float | bool]: times = np.asarray(job.results.get_history_property("Time", "MDHistory"), dtype=float) energies = np.asarray(job.results.get_history_property("PotentialEnergy", "MDHistory"), dtype=float) temperatures = np.asarray(job.results.get_history_property("Temperature", "MDHistory"), dtype=float) start = float(times.min() + 0.5 * (times.max() - times.min())) mask = times >= start trailing_times = times[mask] trailing_energies = energies[mask] trailing_temperatures = temperatures[mask] if len(trailing_times) < 10: raise RuntimeError("Too few trajectory samples for the five-block equilibration test") blocks = np.array_split(np.arange(len(trailing_times)), 5) block_times = np.asarray([trailing_times[block].mean() for block in blocks]) block_energies = np.asarray([trailing_energies[block].mean() for block in blocks]) slope, _ = np.polyfit(block_times, block_energies, 1) fitted_change = float(slope * (block_times[-1] - block_times[0])) trajectory_std = float(np.std(trailing_energies, ddof=1)) energy_ratio = abs(fitted_change) / trajectory_std if trajectory_std > 0.0 else float("inf") mean_temperature = float(np.mean(trailing_temperatures)) temperature_relative_error = abs(mean_temperature - TARGET_TEMPERATURE_K) / TARGET_TEMPERATURE_K energy_ok = bool(energy_ratio <= 0.5) temperature_ok = bool(temperature_relative_error <= 0.02) return { "window_start_fs": start, "window_end_fs": float(times.max()), "fitted_change_hartree": fitted_change, "trajectory_std_hartree": trajectory_std, "energy_change_in_std": energy_ratio, "mean_temperature_k": mean_temperature, "temperature_relative_error": temperature_relative_error, "energy_ok": energy_ok, "temperature_ok": temperature_ok, "passed": bool(energy_ok and temperature_ok), } def log_equilibration_check(duration_ps: int, metrics: Dict[str, float | bool], streak: int) -> None: log(f"Equilibration check at {duration_ps} ps uses the trailing window " f"{float(metrics['window_start_fs']) / 1000:.3f}-{float(metrics['window_end_fs']) / 1000:.3f} ps.") log(f"Potential-energy criterion: |fitted change| / trajectory standard deviation = " f"{float(metrics['energy_change_in_std']):.4f}; required <= 0.5000; " f"{'PASS' if metrics['energy_ok'] else 'FAIL'}.") log(f"Temperature criterion: trailing mean = {float(metrics['mean_temperature_k']):.3f} K; " f"relative error = {100 * float(metrics['temperature_relative_error']):.3f}%; required <= 2.000%; " f"{'PASS' if metrics['temperature_ok'] else 'FAIL'}.") log(f"Combined equilibration check: {'PASS' if metrics['passed'] else 'FAIL'}; " f"consecutive passing checks = {streak}; required = 2.") def restart_settings(previous_job: AMSJob, target_steps: int) -> Settings: settings = previous_job.settings.copy() md = settings.input.ams.MolecularDynamics if "InitialVelocities" in md: del md.InitialVelocities md.NSteps = target_steps md.Restart = previous_job.results.rkfpath() md.CopyRestartTrajectory = "Yes" settings.input.ams.LoadSystem = [ Settings(_h="", file=previous_job.results.rkfpath(), section="Molecule") ] return settings def read_job_text(path: Path) -> str: chunks: List[str] = [] for suffix in ("*.err", "*.out", "*.log"): for file_path in path.glob(suffix): try: chunks.append(file_path.read_text(errors="replace")) except OSError: pass return "\n".join(chunks) def fragmentation_failure_paths() -> List[Path]: from scm.plams import AMSMSDJob failures: List[Path] = [] for workdir in ROOT.glob("03-production_workdir*"): for path in workdir.glob("msd_*ps"): if not path.is_dir(): continue try: if AMSMSDJob.load_external(str(path)).ok(): continue except Exception: pass text = read_job_text(path).lower() if "bond" in text and any(word in text for word in ("change", "changed", "molecule", "fragment")): failures.append(path) return sorted(failures, key=lambda path: natural_key(str(path)))