diff --git a/arc/job/adapters/ase_adapter.py b/arc/job/adapters/ase_adapter.py index 740d3498ea..117f0a69f6 100644 --- a/arc/job/adapters/ase_adapter.py +++ b/arc/job/adapters/ase_adapter.py @@ -202,6 +202,11 @@ def write_input_file(self) -> None: 'irc_direction': self.irc_direction, 'settings': self.determine_settings(), } + if self.job_type == 'scan': + # A 'scan' job sweeps a whole rotor in one process: pass the torsion(s) and the + # scan resolution so ase_script.py can run the relaxed torsional scan itself. + input_dict['torsions'] = self.torsions + input_dict['scan_res'] = self.scan_res save_yaml_file(os.path.join(self.local_path, 'input.yml'), input_dict) def warn_if_unreliable_uma_sp(self) -> bool: diff --git a/arc/job/adapters/ase_test.py b/arc/job/adapters/ase_test.py index d0e6a542fc..9ec1480998 100644 --- a/arc/job/adapters/ase_test.py +++ b/arc/job/adapters/ase_test.py @@ -12,10 +12,30 @@ from unittest.mock import patch import numpy as np +from ase import Atoms +from ase.calculators.emt import EMT + from arc.common import ARC_TESTING_PATH, read_yaml_file, save_yaml_file from arc.job.adapters.ase_adapter import ASEAdapter +from arc.parser.parser import parse_1d_scan_energies from arc.species.species import ARCSpecies -from arc.job.adapters.scripts.ase_script import to_kJmol, numpy_vibrational_analysis, is_linear +from arc.job.adapters.scripts.ase_script import (is_linear, + numpy_vibrational_analysis, + relaxed_torsion_scan, + rotor_top, + run_torsion_scan, + to_kJmol) + +ETHANE_XYZ = {'symbols': ('C', 'C', 'H', 'H', 'H', 'H', 'H', 'H'), + 'isotopes': (12, 12, 1, 1, 1, 1, 1, 1), + 'coords': ((0.0, 0.0, 0.761395), + (0.0, 0.0, -0.761395), + (0.0, 1.017111, 1.157241), + (0.880844, -0.508556, 1.157241), + (-0.880844, -0.508556, 1.157241), + (0.880844, 0.508556, -1.157241), + (-0.880844, 0.508556, -1.157241), + (0.0, -1.017111, -1.157241))} class TestASEAdapter(unittest.TestCase): @@ -204,6 +224,96 @@ def test_numpy_vibrational_analysis(self): for i, val in enumerate(freqs[5:]): self.assertAlmostEqual(results['freqs'][i], val, delta=1e-3) + def test_write_input_file_scan(self): + """Test that a scan job writes torsions and scan resolution into the input file""" + scan_dir = os.path.join(self.project_directory, 'scan_input') + os.makedirs(scan_dir, exist_ok=True) + job = ASEAdapter(execution_type='incore', + job_type='scan', + project='test_scan', + project_directory=scan_dir, + species=[ARCSpecies(label='ethane', xyz=ETHANE_XYZ)], + torsions=[[2, 0, 1, 5]], + args={'keyword': {'calculator': 'uma', 'model': 'uma-s-1p2'}, + 'trsh': {'scan_res': 8.0}}, + testing=True) + job.local_path = scan_dir + job.write_input_file() + data = read_yaml_file(os.path.join(scan_dir, 'input.yml')) + self.assertEqual(data['job_type'], 'scan') + self.assertEqual(data['torsions'], [[2, 0, 1, 5]]) + self.assertEqual(data['scan_res'], 8.0) + + def test_rotor_top(self): + """Test resolving the rotating group across a pivot bond""" + atoms = Atoms(symbols=ETHANE_XYZ['symbols'], positions=ETHANE_XYZ['coords']) + # The dihedral is [2, 0, 1, 5]; pivots are atoms 0 and 1. The group on the atom-1 side is + # that carbon and its three hydrogens (5, 6, 7). + self.assertEqual(rotor_top(atoms, 0, 1), [1, 5, 6, 7]) + self.assertEqual(rotor_top(atoms, 1, 0), [0, 2, 3, 4]) + # A pivot bond inside a ring is ill-defined for a 1D rotor. + cyclopropane = Atoms(symbols=('C', 'C', 'C'), + positions=((0.0, 0.87, 0.0), (0.75, -0.43, 0.0), (-0.75, -0.43, 0.0))) + with self.assertRaises(ValueError): + rotor_top(cyclopropane, 0, 1) + # Pivots that are not bonded cannot define a rotating group. + with self.assertRaises(ValueError): + rotor_top(atoms, 2, 5) + + def test_relaxed_torsion_scan(self): + """Test a full 1D relaxed torsional scan on the machinery (EMT keeps it hermetic)""" + atoms = Atoms(symbols=ETHANE_XYZ['symbols'], positions=ETHANE_XYZ['coords']) + atoms.calc = EMT() + result = relaxed_torsion_scan(atoms, [2, 0, 1, 5], step_deg=8.0, nsteps=45, + fmax=0.05, steps=100) + # 8 deg over 360 deg is 46 points including the duplicated endpoint. + self.assertEqual(len(result['energies']), 46) + self.assertEqual(len(result['angles']), 46) + self.assertEqual(result['angles'][0], 0.0) + self.assertEqual(result['angles'][-1], 360.0) + # The endpoint duplicates the start (a periodic grid). + self.assertEqual(result['energies'][0], result['energies'][-1]) + self.assertTrue(all(np.isfinite(e) for e in result['energies'])) + self.assertEqual(result['top'], [1, 5, 6, 7]) + self.assertIn('fmax_worst', result) + self.assertIn('branch_gap_max', result) + # A non-positive resolution or a zero-length grid is rejected before any division. + with self.assertRaises(ValueError): + relaxed_torsion_scan(atoms, [2, 0, 1, 5], step_deg=8.0, nsteps=0) + with self.assertRaises(ValueError): + relaxed_torsion_scan(atoms, [2, 0, 1, 5], step_deg=0.0, nsteps=45) + + def test_run_torsion_scan_output_is_parseable(self): + """Test that run_torsion_scan output is readable by ARC's YAML scan parser, in kJ/mol""" + atoms = Atoms(symbols=ETHANE_XYZ['symbols'], positions=ETHANE_XYZ['coords']) + atoms.calc = EMT() + input_dict = {'torsions': [[2, 0, 1, 5]], 'scan_res': 8.0} + result = run_torsion_scan(atoms, input_dict, settings={'fmax': 0.05, 'steps': 100}) + self.assertEqual(len(result['energies']), 46) + scan_dir = os.path.join(self.project_directory, 'scan_output') + os.makedirs(scan_dir, exist_ok=True) + out_path = os.path.join(scan_dir, 'output.yml') + save_yaml_file(out_path, {'energies': result['energies'], 'angles': result['angles']}) + energies, angles = parse_1d_scan_energies(log_file_path=out_path) + self.assertIsNotNone(energies) + self.assertIsNotNone(angles) + self.assertEqual(len(energies), 46) + self.assertEqual(len(angles), 46) + self.assertAlmostEqual(min(energies), 0.0) # parser zeroes the minimum, in kJ/mol + self.assertGreaterEqual(min(energies), 0.0) + + def test_run_torsion_scan_rejects_multiple_torsions(self): + """Test that a 1D scan refuses more than one torsion""" + atoms = Atoms(symbols=ETHANE_XYZ['symbols'], positions=ETHANE_XYZ['coords']) + atoms.calc = EMT() + with self.assertRaises(ValueError): + run_torsion_scan(atoms, {'torsions': [[2, 0, 1, 5], [3, 0, 1, 6]], 'scan_res': 8.0}, settings={}) + with self.assertRaises(ValueError): + run_torsion_scan(atoms, {'torsions': None, 'scan_res': 8.0}, settings={}) + # A non-positive scan resolution is rejected before the 360/scan_res division. + with self.assertRaises(ValueError): + run_torsion_scan(atoms, {'torsions': [[2, 0, 1, 5]], 'scan_res': 0.0}, settings={}) + @classmethod def tearDownClass(cls): """ diff --git a/arc/job/adapters/scripts/ase_script.py b/arc/job/adapters/scripts/ase_script.py index 2612b94de9..f0b1905b32 100644 --- a/arc/job/adapters/scripts/ase_script.py +++ b/arc/job/adapters/scripts/ase_script.py @@ -15,6 +15,7 @@ from ase import Atoms from ase.constraints import FixInternals +from ase.neighborlist import build_neighbor_list, natural_cutoffs from ase.optimize import BFGS, LBFGS, GPMin from ase.optimize.sciopt import SciPyFminBFGS, SciPyFminCG from ase.vibrations import Vibrations @@ -130,6 +131,205 @@ def apply_constraints(atoms: Atoms, constraints_data: list): dihedrals.append([constraint[1], indices]) atoms.set_constraint(FixInternals(bonds=bonds, angles_deg=angles, dihedrals_deg=dihedrals)) + +def rotor_top(atoms: Atoms, pivot_1: int, pivot_2: int) -> list: + """ + Determine the rotating group (the "top") on the ``pivot_2`` side of the ``pivot_1``-``pivot_2`` bond. + + The two pivots are the central atoms of the scanned torsion; breaking their bond splits the + molecule into two fragments, and this returns the atom indices reachable from ``pivot_2``. + + This is the same "top" that ``arc.common.determine_top_group_indices`` computes, but kept as a + standalone reimplementation on purpose: this script runs in the calculator's own environment + (e.g. ``uma_env``), which has neither ARC nor RMG-Py, so it can import neither that function nor + the RMG ``Molecule`` it operates on. Connectivity is instead rebuilt from the geometry via ASE's + neighbor list. Keep the two traversals in agreement if either changes. + + Args: + atoms (Atoms): The molecule. + pivot_1 (int): The 0-indexed atom on the fixed side of the rotation axis. + pivot_2 (int): The 0-indexed atom on the rotating side of the rotation axis. + + Raises: + ValueError: If the pivots are not bonded, or if the pivot bond is part of a ring (a 1D + rotor is then ill-defined). + + Returns: + list: The sorted 0-indexed atoms of the rotating group. + """ + nl = build_neighbor_list(atoms, natural_cutoffs(atoms, mult=1.2), + self_interaction=False, bothways=True) + adj = {i: set(nl.get_neighbors(i)[0]) for i in range(len(atoms))} + if pivot_2 not in adj[pivot_1]: + raise ValueError(f'The pivot atoms {pivot_1} and {pivot_2} are not bonded; ' + f'cannot determine the rotating group.') + adj[pivot_1].discard(pivot_2) + adj[pivot_2].discard(pivot_1) + seen, stack = {pivot_2}, [pivot_2] + while stack: + cur = stack.pop() + for nb in adj[cur]: + if nb not in seen: + seen.add(nb) + stack.append(nb) + if pivot_1 in seen: + raise ValueError(f'The pivot bond {pivot_1}-{pivot_2} is part of a ring; ' + f'a 1D rotor scan is ill-defined.') + return sorted(int(x) for x in seen) + + +def _scan_walk(atoms: Atoms, torsion: list, step_deg: float, nsteps: int, top: list, + fmax: float, steps: int, direction: int, optimizer) -> tuple: + """ + Perform one directional sequential relaxed torsional walk. + + Each grid point starts from the previous point's relaxed geometry (a hysteretic walk), so a + coupled coordinate can carry the walk into a different conformer and never return - which is + why a full scan runs this in both directions and keeps the lower branch at each point. + + Args: + atoms (Atoms): The molecule with a calculator attached (not modified; a copy is walked). + torsion (list): The 0-indexed four atoms (i, j, k, l) defining the scanned dihedral. + step_deg (float): The dihedral increment in degrees. + nsteps (int): The number of increments (the walk visits ``nsteps`` + 1 points). + top (list): The 0-indexed rotating group. + fmax (float): The force convergence criterion in eV/Angstrom. + steps (int): The maximum optimizer steps per point. + direction (int): +1 for a 0->360 walk, -1 for a 360->0 walk. + optimizer: The ASE optimizer class to relax each point. + + Returns: + tuple: (energies in eV, residual max-forces in eV/Angstrom), each of length ``nsteps`` + 1. + """ + i, j, k, l = torsion + work = atoms.copy() + work.calc = atoms.calc + energies, residual_forces = list(), list() + for n in range(nsteps + 1): + if n: + work.rotate_dihedral(i, j, k, l, angle=direction * step_deg, indices=top) + target = work.get_dihedral(i, j, k, l) + work.set_constraint(FixInternals(dihedrals_deg=[[target, [i, j, k, l]]])) + opt = optimizer(work, logfile=None) + opt.run(fmax=fmax, steps=steps) + # Residual force WITH the dihedral constraint still applied is the convergence criterion: + # reading opt.converged() after clearing the constraint re-evaluates the free torsional + # force, which is nonzero everywhere except the stationary points. + residual = float(np.sqrt((work.get_forces() ** 2).sum(axis=1)).max()) + work.set_constraint() + energies.append(float(work.get_potential_energy())) + residual_forces.append(residual) + return energies, residual_forces + + +def relaxed_torsion_scan(atoms: Atoms, torsion: list, step_deg: float, nsteps: int, + top: list = None, fmax: float = 0.005, steps: int = 500, + optimizer=LBFGS) -> dict: + """ + Perform a full 1D relaxed torsional scan in a single process. + + The torsion is swept 0->360 and 360->0 from the same starting geometry, and the lower energy + of the two walks is kept at each grid point. A single directional walk is path-dependent: if a + coupled coordinate flips partway round, that walk finishes in a different conformer and its + curve is not a torsional potential. Walking both ways and taking the pointwise minimum recovers + the lowest torsional path, which is the one the downstream statistical mechanics wants. + + Points that do not fully converge are kept (never dropped), so the periodic grid the downstream + Fourier fit reads is always complete; ``fmax_worst`` reports the largest residual force. + + Args: + atoms (Atoms): The molecule with a calculator attached. + torsion (list): The 0-indexed four atoms (i, j, k, l) defining the scanned dihedral. + step_deg (float): The dihedral increment in degrees. + nsteps (int): The number of increments; the returned grid has ``nsteps`` + 1 points. + top (list, optional): The 0-indexed rotating group; determined from the graph if not given. + fmax (float): The force convergence criterion in eV/Angstrom. + steps (int): The maximum optimizer steps per point. + optimizer: The ASE optimizer class to relax each point. + + Raises: + ValueError: If ``nsteps`` is less than 1 or ``step_deg`` is not positive. + + Returns: + dict: ``energies`` (absolute, in eV, ``nsteps`` + 1 points with the endpoint duplicating the + start), ``angles`` (degrees, 0..360), ``top``, and convergence diagnostics. + """ + if nsteps < 1: + raise ValueError(f'A relaxed torsion scan needs at least one increment, but got nsteps={nsteps}.') + if step_deg <= 0: + raise ValueError(f'The dihedral increment must be positive, but got step_deg={step_deg}.') + i, j, k, l = torsion + if top is None: + top = rotor_top(atoms, j, k) + + e_f, r_f = _scan_walk(atoms, torsion, step_deg, nsteps, top, fmax, steps, +1, optimizer) + e_b, r_b = _scan_walk(atoms, torsion, step_deg, nsteps, top, fmax, steps, -1, optimizer) + + npts = nsteps # unique grid points; the walk's index nsteps maps back onto index 0 + grid_f = np.full(npts, np.inf) + grid_b = np.full(npts, np.inf) + for n in range(nsteps + 1): + p = n % npts # forward: phi = +n*step + grid_f[p] = min(grid_f[p], e_f[n]) + q = (-n) % npts # backward: phi = -n*step + grid_b[q] = min(grid_b[q], e_b[n]) + merged = np.minimum(grid_f, grid_b) + energies = merged.tolist() + energies.append(energies[0]) # duplicate endpoint: 46-point protocol + angles = [n * step_deg for n in range(nsteps + 1)] + return { + 'energies': energies, + 'angles': angles, + 'top': top, + 'branch_gap_max': float(np.abs(grid_f - grid_b).max()), + 'fmax_worst': float(max(r_f + r_b)), + 'converged': bool(max(r_f + r_b) <= fmax), + } + + +def run_torsion_scan(atoms: Atoms, input_dict: dict, settings: dict) -> dict: + """ + Run a 1D relaxed torsional scan for an ARC ``scan`` job and shape it for ARC's YAML parser. + + Args: + atoms (Atoms): The molecule with a calculator attached. + input_dict (dict): The job input; ``torsions`` (a list holding one 0-indexed torsion) and + ``scan_res`` (the dihedral increment in degrees) are read here. + settings (dict): ASE run settings; optional ``fmax``, ``steps``, ``optimizer`` override the + relaxed-scan defaults (LBFGS, fmax 0.005 eV/Angstrom, 500 steps per point). + + Raises: + ValueError: If no torsion is supplied or more than one is (a 1D scan sweeps a single + torsion), or if ``scan_res`` is not positive. + + Returns: + dict: ``energies`` (Hartree) and ``angles`` (degrees) for ARC's YAML scan parser, plus the + rotating ``top`` and convergence diagnostics. + """ + torsions = input_dict.get('torsions') + if not torsions: + raise ValueError("A 'scan' job requires a torsion, but none was supplied.") + if len(torsions) > 1: + raise ValueError(f"A 1D 'scan' job sweeps a single torsion, but got {len(torsions)}: {torsions}.") + torsion = list(torsions[0]) + scan_res = float(input_dict.get('scan_res', 8.0)) + if scan_res <= 0: + raise ValueError(f"The scan resolution must be positive, but got scan_res={scan_res}.") + nsteps = int(round(360.0 / scan_res)) + fmax = float(settings.get('fmax', 0.005)) + steps = int(settings.get('steps', 500)) + engine_dict = {'bfgs': BFGS, 'lbfgs': LBFGS, 'gpmin': GPMin, + 'scipyfminbfgs': SciPyFminBFGS, 'scipyfmincg': SciPyFminCG} + optimizer = engine_dict.get(str(settings.get('optimizer', 'LBFGS')).lower(), LBFGS) + + result = relaxed_torsion_scan(atoms, torsion, scan_res, nsteps, top=None, + fmax=fmax, steps=steps, optimizer=optimizer) + # ARC's YAML scan parser reads absolute energies and subtracts the minimum before converting + # Hartree->kJ/mol, so the energies must be in Hartree. + result['energies'] = [energy * e / E_h for energy in result['energies']] + return result + + def is_linear(atoms: Atoms) -> bool: """ Determine whether an Atoms object represents a linear molecule. @@ -345,6 +545,13 @@ def save_current_geometry(out_dict, atoms_obj, input_xyz): if job_type == 'sp': output['sp'] = to_kJmol(atoms.get_potential_energy()) + if job_type == 'scan': + try: + output.update(run_torsion_scan(atoms, input_dict, settings)) + except Exception as exc: + output['success'] = False + output['error'] = f"Torsion scan failed: {exc}" + if job_type in ['opt', 'conf_opt', 'optfreq', 'directed_scan']: fmax = float(settings.get('fmax', 0.001)) steps = int(settings.get('steps', 1000))