Add a one-process relaxed rotor scan to the ASE adapter

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:

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):
         """

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))