#!/usr/bin/env amspython from __future__ import annotations from dataclasses import dataclass import numpy as np import scm.plams as plams from scm.base import ChemicalSystem from scm.input_classes import AMS from scm.plams import AMSJob, Settings SMILES = "C=CC(CCCC)C1C=CC(C)=CC1=O" TARGET_DISTANCE_ANGSTROM = 1.5 SCAN_STEP_ANGSTROM = 0.2 @dataclass(frozen=True) class AtomPair: carbonyl_oxygen: int terminal_vinyl_carbon: int def ensure_bonds(system: object) -> None: has_bonds = getattr(system, "has_bonds", None) if callable(has_bonds): if not has_bonds(): system.guess_bonds() return bonds = getattr(system, "bonds", None) if bonds is not None and len(bonds) == 0: system.guess_bonds() def atom_symbol(system: ChemicalSystem, index: int) -> str: return system.atoms[index].element.symbol def adjacency(system: ChemicalSystem) -> dict[int, list[tuple[int, float]]]: graph: dict[int, list[tuple[int, float]]] = {i: [] for i in range(len(system.atoms))} for i, j, bond in system.bonds: graph[i].append((j, bond.order)) graph[j].append((i, bond.order)) return graph def heavy_neighbors(system: ChemicalSystem, graph: dict[int, list[tuple[int, float]]], index: int) -> list[tuple[int, float]]: return [(nbr, order) for nbr, order in graph[index] if atom_symbol(system, nbr) != "H"] def distance_angstrom(system: ChemicalSystem, i: int, j: int) -> float: coords_i = np.array(system.atoms[i].coords) coords_j = np.array(system.atoms[j].coords) return float(np.linalg.norm(coords_i - coords_j)) def find_target_atoms(system: ChemicalSystem) -> AtomPair: ensure_bonds(system) graph = adjacency(system) carbonyl_oxygen_candidates: list[int] = [] for index in range(len(system.atoms)): if atom_symbol(system, index) != "O": continue nbrs = heavy_neighbors(system, graph, index) if len(nbrs) != 1: continue nbr, order = nbrs[0] if atom_symbol(system, nbr) == "C" and order >= 1.5: carbonyl_oxygen_candidates.append(index) if len(carbonyl_oxygen_candidates) != 1: raise RuntimeError(f"Expected exactly one carbonyl oxygen, found {carbonyl_oxygen_candidates}") carbonyl_oxygen = carbonyl_oxygen_candidates[0] terminal_vinyl_candidates: list[int] = [] for index in range(len(system.atoms)): if atom_symbol(system, index) != "C": continue nbrs = heavy_neighbors(system, graph, index) if len(nbrs) != 1: continue nbr, order = nbrs[0] if atom_symbol(system, nbr) == "C" and order >= 1.5: terminal_vinyl_candidates.append(index) if not terminal_vinyl_candidates: raise RuntimeError("Could not find a terminal vinyl carbon from the bond graph") terminal_vinyl_carbon = max( terminal_vinyl_candidates, key=lambda idx: distance_angstrom(system, carbonyl_oxygen, idx), ) return AtomPair(carbonyl_oxygen=carbonyl_oxygen, terminal_vinyl_carbon=terminal_vinyl_carbon) def validate(settings: Settings) -> None: AMS.from_settings(settings) def dftb_settings() -> Settings: settings = Settings() settings.input.DFTB.Model = "GFN1-xTB" return settings def reactant_settings() -> Settings: settings = dftb_settings() settings.input.ams.Task = "GeometryOptimization" settings.input.ams.Properties.NormalModes = "Yes" validate(settings) return settings def pesscan_settings(atom_pair: AtomPair, start_distance: float, end_distance: float, npoints: int) -> Settings: settings = dftb_settings() settings.input.ams.Task = "PESScan" settings.input.ams.PESScan.ScanCoordinate.nPoints = npoints settings.input.ams.PESScan.ScanCoordinate.Distance = ( f"{atom_pair.carbonyl_oxygen + 1:d} {atom_pair.terminal_vinyl_carbon + 1:d} " f"{start_distance:.6f} {end_distance:.6f}" ) validate(settings) return settings def ts_settings() -> Settings: settings = dftb_settings() settings.input.ams.Task = "TransitionStateSearch" settings.input.ams.Properties.NormalModes = "Yes" settings.input.ams.GeometryOptimization.InitialHessian.Type = "Calculate" validate(settings) return settings def product_settings() -> Settings: settings = dftb_settings() settings.input.ams.Task = "GeometryOptimization" settings.input.ams.Properties.NormalModes = "Yes" validate(settings) return settings def run_job(name: str, molecule: ChemicalSystem | dict[str, ChemicalSystem], settings: Settings) -> AMSJob: job = AMSJob(molecule=molecule, settings=settings, name=name) result = job.run() if not result.ok(): raise RuntimeError(f"{name} failed") return job def main() -> None: plams.init(folder="01-run_workdir") system = ChemicalSystem.from_smiles(SMILES) ensure_bonds(system) initial_atom_pair = find_target_atoms(system) print( "Initial atom indices (1-based):", f"carbonyl O={initial_atom_pair.carbonyl_oxygen + 1}", f"terminal vinyl C={initial_atom_pair.terminal_vinyl_carbon + 1}", ) reactant_job = run_job("reactant_opt", system, reactant_settings()) reactant_system = reactant_job.results.get_main_system() ensure_bonds(reactant_system) atom_pair = find_target_atoms(reactant_system) start_distance = distance_angstrom(reactant_system, atom_pair.carbonyl_oxygen, atom_pair.terminal_vinyl_carbon) if start_distance <= TARGET_DISTANCE_ANGSTROM: raise RuntimeError( f"Reactant O...C distance is already {start_distance:.3f} A, cannot scan down to {TARGET_DISTANCE_ANGSTROM:.3f} A" ) npoints = int(round((start_distance - TARGET_DISTANCE_ANGSTROM) / SCAN_STEP_ANGSTROM)) + 1 if npoints < 2: raise RuntimeError(f"Invalid PES scan point count from start distance {start_distance:.3f} A") pesscan_job = run_job( "pes_scan", reactant_system, pesscan_settings(atom_pair, start_distance, TARGET_DISTANCE_ANGSTROM, npoints), ) pesscan_results = pesscan_job.results.get_pesscan_results() energies = np.array(pesscan_results["PES"], dtype=float) distances = np.array(pesscan_results["RaveledPESCoords"][0], dtype=float) molecules = pesscan_results["Molecules"] highest_index = int(np.argmax(energies)) shortest_index = int(np.argmin(distances)) highest_system = molecules[highest_index] product_guess = molecules[shortest_index] ensure_bonds(highest_system) ensure_bonds(product_guess) ts_job = run_job( "ts_search", {"": highest_system, "final": product_guess}, ts_settings(), ) ts_frequencies = np.array(ts_job.results.get_frequencies(unit="cm^-1"), dtype=float) n_imaginary = int(np.count_nonzero(ts_frequencies < 0.0)) print(f"Transition state imaginary frequencies: {n_imaginary}") if n_imaginary != 1: raise RuntimeError(f"Transition state verification failed: expected 1 imaginary frequency, found {n_imaginary}") product_job = run_job("product_opt", product_guess, product_settings()) product_system = product_job.results.get_main_system() ensure_bonds(product_system) print("Finished workflow") print( f"Final atom indices (1-based) used for scan: O={atom_pair.carbonyl_oxygen + 1}, " f"C={atom_pair.terminal_vinyl_carbon + 1}" ) print(f"Scan start distance: {start_distance:.3f} A") print(f"Requested scan end distance: {TARGET_DISTANCE_ANGSTROM:.3f} A") print(f"Number of PES points: {npoints}") if __name__ == "__main__": main()