Skip to content

CREST Adapter#807

Open
calvinp0 wants to merge 241 commits into
mainfrom
crest_adapter
Open

CREST Adapter#807
calvinp0 wants to merge 241 commits into
mainfrom
crest_adapter

Conversation

@calvinp0

@calvinp0 calvinp0 commented Nov 27, 2025

Copy link
Copy Markdown
Member

Addition of CREST Adapter that complements the heuristic adapter.

This pull request adds support for the CREST conformer and transition state (TS) search method to the ARC project, along with several related improvements and code cleanups. The most important changes include integrating CREST as a TS search adapter, updating configuration and constants, and enhancing the heuristics TS search logic for better provenance tracking and code clarity.

CREST Integration:

  • Added CREST as a supported TS search method: updated JobEnum (arc/job/adapter.py), included CREST in the list of adapters and RMG family mapping, and registered it as a default incore adapter (arc/job/adapters/common.py, arc/job/adapters/ts/__init__.py). [1] [2] [3] [4]
  • Implemented a new test suite for CREST input generation (arc/job/adapters/ts/crest_test.py).
  • Added a Makefile target and installation script for CREST (Makefile). [1] [2]

Constants and Configuration:

  • Added the angstrom_to_bohr conversion constant to both Cython and Python constants modules (arc/constants.pxd, arc/constants.py). [1] [2] [3]

Heuristics TS Search Enhancements and Refactoring:

  • Refactored heuristics TS search logic to track and combine method provenance for TS guesses, allowing for more precise attribution when multiple methods contribute to a guess (arc/job/adapters/ts/heuristics.py). [1] [2] [3] [4]
  • Improved code readability and maintainability by reformatting imports and function calls, and clarifying data structures and comments in heuristics TS search (arc/job/adapters/ts/heuristics.py). [1] [2] [3] [4] [5] [6] [7]
    .

Comment thread arc/job/adapters/ts/crest.py Fixed
Comment thread arc/job/adapters/ts/crest.py Fixed
Comment thread arc/job/adapters/ts/crest.py Fixed
Comment thread arc/settings/settings.py Fixed
Comment thread arc/settings/settings.py Fixed
Comment thread arc/job/adapters/ts/autotst_ts.py Fixed
Comment thread arc/job/adapters/ts/crest.py Fixed
Comment thread arc/job/adapters/ts/heuristics.py Fixed
Comment thread arc/job/adapters/ts/heuristics.py Fixed
Comment thread arc/job/adapters/ts/heuristics.py Fixed
Comment thread arc/job/adapters/ts/crest.py Fixed

This comment was marked as resolved.

Comment thread arc/job/adapters/ts/crest.py Fixed

This comment was marked as resolved.

Comment thread arc/job/adapters/ts/heuristics.py Fixed
Comment thread arc/settings/settings.py Fixed
Comment thread arc/settings/settings.py Fixed
@calvinp0
calvinp0 force-pushed the crest_adapter branch 2 times, most recently from 3e88c36 to 7674f5f Compare February 2, 2026 09:49
@calvinp0
calvinp0 requested a review from Copilot February 2, 2026 12:10

This comment was marked as resolved.

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copilot encountered an error and was unable to review this pull request. You can try again by re-requesting a review.

This comment was marked as resolved.

@calvinp0
calvinp0 force-pushed the crest_adapter branch 2 times, most recently from 9be935e to eab7647 Compare February 4, 2026 20:14
@calvinp0
calvinp0 requested a review from alongd February 4, 2026 21:30
Register 'Disproportionation' in both TS-adapter gates for autotst:
the adapter's own supported_families and ts_adapters_by_rmg_family.
intra_H_migration was already enabled and is untouched.

Note: Disproportionation requires the AutoTST env on a branch whose
SUPPORTED_FAMILIES includes it (e.g. fix/disproportionation-support),
not AutoTST main. Adds a focused test asserting both gates admit it.
shutil.rmtree(path, ignore_errors=True)
try:
os.rmdir(self.project_dir)
except OSError:

Check notice

Code scanning / CodeQL

Empty except Note

'except' clause does nothing but pass and there is no explanatory comment.

Copilot Autofix

AI 8 days ago

The best fix is to replace the empty except with explicit handling that preserves current functionality (best-effort cleanup) while making intent clear and avoiding silent swallowing.

In arc/job/adapters/ts/autotst_ts_test.py, update _remove_test_dir so that:

  • os.rmdir(self.project_dir) remains in a try.
  • except OSError becomes except OSError as e.
  • Add an explanatory comment and only suppress expected race/cleanup conditions (directory not empty / not found), re-raising unexpected OS errors.

This keeps tests robust under parallel execution and satisfies the static-analysis concern by avoiding an empty handler.

Suggested changeset 1
arc/job/adapters/ts/autotst_ts_test.py

Autofix patch

Autofix patch
Run the following command in your local git repository to apply this patch
cat << 'EOF' | git apply
diff --git a/arc/job/adapters/ts/autotst_ts_test.py b/arc/job/adapters/ts/autotst_ts_test.py
--- a/arc/job/adapters/ts/autotst_ts_test.py
+++ b/arc/job/adapters/ts/autotst_ts_test.py
@@ -44,8 +44,12 @@
         shutil.rmtree(path, ignore_errors=True)
         try:
             os.rmdir(self.project_dir)
-        except OSError:
-            pass
+        except OSError as e:
+            # Best-effort cleanup: under parallel execution the shared parent directory
+            # may already be removed or still contain other tests' subdirectories.
+            if not os.path.isdir(self.project_dir):
+                return
+            raise e
 
     def get_adapter(self, dir_name: str) -> AutoTSTAdapter:
         """A helper function to instantiate an AutoTSTAdapter instance."""
EOF
@@ -44,8 +44,12 @@
shutil.rmtree(path, ignore_errors=True)
try:
os.rmdir(self.project_dir)
except OSError:
pass
except OSError as e:
# Best-effort cleanup: under parallel execution the shared parent directory
# may already be removed or still contain other tests' subdirectories.
if not os.path.isdir(self.project_dir):
return
raise e

def get_adapter(self, dir_name: str) -> AutoTSTAdapter:
"""A helper function to instantiate an AutoTSTAdapter instance."""
Copilot is powered by AI and may make mistakes. Always verify output.
calvinp0 added 7 commits July 15, 2026 12:51
An impossible spin multiplicity (e.g. an even multiplicity for a species
with an even electron count) was forwarded all the way to the ESS, where
it surfaced as a cryptic non-retryable Gaussian GL301 InputError and
aborted the species after wasting a full run.

Enforce the parity relation
    (total_electrons - net_charge) % 2 == (multiplicity - 1) % 2
in ARCSpecies once charge and multiplicity are finalized, raising a clear
SpeciesError that names the label, electron count, net charge, requested
multiplicity, and the nearest valid multiplicities. This is an inviolable
parity relation, so the guard has zero false positives. It is skipped
gracefully when the electron count is not yet known (e.g. a TS or an
xyz-less species with no perceived Molecule).
A fine transition-state optimization that dead-ends on the Gaussian l9999
"Optimization stopped" / "Number of steps exceeded" oscillation (classified
as ['MaxOptCycles', 'GL9999']) was previously escalated only through the
Hessian/step ladder (maxcycle=200 -> recalcfc=5 -> calcall -> RFO) and then
abandoned via 'all_attempted'. opt=(cartesian) - wired until now only to the
GL103/InternalCoordinateError path - was never tried, even though the failure
is a redundant-internal-coordinate oscillation that Cartesian coordinates
sidestep. This lost chemically-sensible TS guesses (e.g. arcbench reaction_16
C=C + HCl, heuristics-xy conformers 0/1).

Insert opt=(cartesian) as an early rung in the MaxOptCycles ladder, right
after the cheap maxcycle bump and before the expensive Hessian/algorithm
escalation, gated tightly to the fine TS opt (is_ts and fine). Ground-state
opts (is_ts=False) and coarse TS opts (fine=False) are unchanged, as is the
GL103 path. Cartesian is recorded once (bare 'cartesian' marker, consistent
with trsh_keyword_cartesian) and folded into the merged opt route so it is
carried through the remaining rungs; the ladder still terminates via
'all_attempted' if Cartesian also fails.

Also fixes a latent cartesian-loss in mixed GL103+MaxOptCycles histories,
where the opt-route rewrite would previously clobber a standalone
opt=(cartesian).
seed = {'xyz': xyz, 'family': rxn.family}
constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed)
self.assertIsInstance(constraints, dict)
self.assertTrue({'A', 'H', 'B', 'atoms', 'distance_pairs', 'angle_atoms'} <= set(constraints))

Check notice

Code scanning / CodeQL

Imprecise assert Note

assertTrue(a <= b) cannot provide an informative message. Using assertLessEqual(a, b) instead will give more informative messages.

Copilot Autofix

AI 4 days ago

Replace the imprecise boolean assertion with a specific unittest assertion tailored to subset comparison.

Best fix in this file: in arc/job/adapters/ts/heuristics_test.py, method test_get_wrapper_constraints_crest, replace:

  • self.assertTrue(required_keys <= set(constraints)) style check
    with:
  • self.assertLessEqual(required_keys, set(constraints))

This preserves exact test logic (subset relationship) while improving failure output (you’ll see the compared sets). No new imports or helper methods are needed.

Suggested changeset 1
arc/job/adapters/ts/heuristics_test.py

Autofix patch

Autofix patch
Run the following command in your local git repository to apply this patch
cat << 'EOF' | git apply
diff --git a/arc/job/adapters/ts/heuristics_test.py b/arc/job/adapters/ts/heuristics_test.py
--- a/arc/job/adapters/ts/heuristics_test.py
+++ b/arc/job/adapters/ts/heuristics_test.py
@@ -2391,7 +2391,7 @@
         seed = {'xyz': xyz, 'family': rxn.family}
         constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed)
         self.assertIsInstance(constraints, dict)
-        self.assertTrue({'A', 'H', 'B', 'atoms', 'distance_pairs', 'angle_atoms'} <= set(constraints))
+        self.assertLessEqual({'A', 'H', 'B', 'atoms', 'distance_pairs', 'angle_atoms'}, set(constraints))
         self.assertEqual(
             (constraints['A'], constraints['H'], constraints['B']),
             constraints['angle_atoms'],
EOF
@@ -2391,7 +2391,7 @@
seed = {'xyz': xyz, 'family': rxn.family}
constraints = get_wrapper_constraints(wrapper='crest', reaction=rxn, seed=seed)
self.assertIsInstance(constraints, dict)
self.assertTrue({'A', 'H', 'B', 'atoms', 'distance_pairs', 'angle_atoms'} <= set(constraints))
self.assertLessEqual({'A', 'H', 'B', 'atoms', 'distance_pairs', 'angle_atoms'}, set(constraints))
self.assertEqual(
(constraints['A'], constraints['H'], constraints['B']),
constraints['angle_atoms'],
Copilot is powered by AI and may make mistakes. Always verify output.
calvinp0 added 12 commits July 15, 2026 19:41
…ny cores

On small double-radical PESs (e.g. OH + OH -> H2O + O) the CREST H-abstraction
$constrain block pinned only the A-H and H-B distances plus a weak positional
restraint, leaving the A...B heavy-heavy separation and the A-H-B angle free. GFN2-xTB
metadynamics then collapsed the linear A...H...B seed into a bent A-B minimum, which the
post-CREST angle guard correctly rejected.

Part 1: for the three-center H-abstraction path, also emit 'distance: A, B, auto' and
'angle: A, H, B, auto' in the $constrain block (xtb supports the angle keyword; see
xtb_adapter.py). The XY four-center path (no angle_atoms) is left unchanged.

Part 2: skip the CREST guess when the reactive core spans essentially the whole molecule
(H_Abstraction with <=1 spectator atom, i.e. <=4 atoms total), where CREST has no
conformational DOF to sample; defer to the other TS-search methods.
An l202 'OptOrientation' error (standard orientation / point group changed
mid-optimization) is a symmetry glitch that nosymm cures - the same remedy
l101/l103 already use. It was mis-grouped with genuine dead-ends in
GAUSSIAN_NON_RETRYABLE_KEYWORDS, which short-circuited trsh_ess_job into a
non-retryable refusal before the Gaussian remedy branch was reached.

Remove 'OptOrientation' from GAUSSIAN_NON_RETRYABLE_KEYWORDS and add 'NoSymm'
to the l202 keyword list in determine_ess_status, so trsh_keyword_nosymm fires
and injects nosymm exactly once. If nosymm also fails the job terminates via
all_attempted (no infinite retry). l101/l103 and other non-retryable keywords
are unchanged.
For a reaction whose family was discovered in the reverse direction, the
first-orientation map_rxn attempt fails on get_template_product_order (or on a
missing template map) and returns None, after which map_reaction routinely
recovers via its existing flip fallback and returns a valid map. The two
intermediate logger.error calls in map_rxn therefore fired as false alarms on
reactions that mapped successfully.

Downgrade both to logger.debug. Control flow is unchanged (both branches still
return None as before), so the returned atom map is identical for every
reaction. A genuine total-mapping failure (both orientations exhausted) is
still surfaced at error level by ARCReaction.atom_map ("could not be atom
mapped").

Add a focused test that drives the first-orientation failure via map_rxn and
asserts no ERROR is logged, while map_reaction still returns a valid map.
ARCSpecies.set_dihedral() correctly refuses to rotate a dihedral that
spans a linear segment (e.g. a collinear O=C=C cumulene) and returns
None. This is a benign no-op, but it was logged at WARNING level and
the mapping engine's backbone-alignment loop re-attempted the same
linear torsion on every iteration, producing repeated warnings.

- set_dihedral(): downgrade the linear-segment message to logger.debug.
  The is_angle_linear guard and the return None are unchanged.
- get_backbone_dihedral_angles(): skip backbone torsions that span a
  linear segment (new is_torsion_linear helper, mirroring the guard) so
  they are never fed into the set_dihedral alignment loop. This is a
  no-op for non-linear torsions and does not change alignment results.

Adds focused tests for both behaviors.
Molpro's `memory,Total=N,m` card is allocated PER MPI PROCESS and the job
runs `molpro -n {cpu_cores}`, but set_input_file_memory hard-coded N as
job_memory_gb * 31.25 (== 125/4, an assumed 4 processes) and never divided
by the actual cpu_cores. At higher rank counts the card demanded several
times the granted node memory -> 'Insufficient memory to allocate'.

(a) Compute the card per-process: ceil(job_memory_gb * 125 / cpu_cores),
    mirroring Orca/CFOUR/TeraChem. 125 MW == 1 GB, so total memory equals
    job_memory_gb for any rank count (previously only correct at 4 ranks).

(b)/(c) Molpro 'Memory' trsh now bounds the requested total by the node
    cap (node_mem * job_max_server_node_memory_allocation) instead of
    inflating unboundedly, and once at the cap reduces the MPI rank count
    (halving cpu_cores, already returned/threaded to set_cpu_and_mem) so the
    per-process share grows without exceeding node memory. When it cannot
    reduce ranks further it emits a terminal error rather than resubmitting
    identically until max_ess_trsh.

Update molpro_test set_input_file_memory expectations (48->37, 8->219,
1->1750 for job_memory_gb=14) and the written memory cards (438->37 at 48
ranks); update/extend trsh_test for the bounded + rank-reducing branch.
The at-cap rank-halving branch set memory = min(memory_gb, max_mem_allocation),
which kept the pre-cap total whenever memory_gb was still below the cap. A job
that needs ~the full node (e.g. 62 GB, cap 60.8) was then declared dead at a
single rank having never actually been offered the cap. Pin the total to the
node cap (memory = max_mem_allocation) while reducing ranks so the per-rank
share is raised against the full node allocation before giving up.

Update trsh_test expectations: the at-cap halving cases now expect memory ==
max_mem_allocation (the cap) at each step; the focused test's at-cap case uses
memory_gb below the cap so the pin lifts the total up to the cap.
…ount)

Atomic nitrogen has 7 electrons, so its multiplicity must be even.
Ground-state N is 4S (multiplicity=4). The old multiplicity=3 is
physically impossible and now raises SpeciesError under the
electron-parity guard, erroring the whole molpro_test class.
Also update the job_7 expected Molpro input wf,spin cards from
spin=2 to spin=3 (Molpro spin = multiplicity - 1).
…d trsh

The prior fix divided the memory,Total=N,m card by cpu_cores on the
assumption it is per-process. molpro26 on zeus proves it is NODE-TOTAL:
a written Total=438 printed 'Total memory per node: 438 MW' while giving
only 'Memory per process: 30 MW' (Molpro reserves ~25% for Global Arrays
and splits the rest across the -n ranks itself). Dividing by cpu_cores
double-divides and starves the node; a 42 GB/12-rank job got 3.5 GB.

- molpro.py set_input_file_memory: N = ceil(job_memory_gb * 125), so the
  node-total card equals the PBS reservation. No cpu_cores division.
- trsh.determine_ess_status: the 'A further ... Increase memory to M'
  branch now captures the target M (last number) as a per-process floor;
  fixes the always-false 'error not in locals()' guard (-> 'not error').
- trsh.trsh_ess_job molpro Memory branch: when a per-process target is
  parsed, size the node-total card as ceil(M*1.5*nprocs/0.822)/125 GB so
  each rank gets >= M MW; keep the memory_gb*3 fallback and the at-cap
  rank-halving. Fixed stray /128 -> node-total sizing.
- Tests: node-total card == ceil(job_memory_gb*125), independent of
  cpu_cores; added a targeted-retry test for the per-process-target path.
The heuristic TS seed-builder (combine_coordinates_with_redundant_atoms)
failed to produce any guess when the H-abstraction acceptor has a linear /
cumulene reactive center (e.g. HCCO + THF <=> CH2CO + THF_rad, where the
ketene product has a collinear O=C=C backbone).

Root cause: in get_new_zmat_2_map, after remove_zmat_atom_0 the reduced
zmat_2 of a linear/cumulene fragment carries a collinear (~0 deg) bond angle.
zmat_to_xyz cannot place a dihedral that references three collinear atoms
(cross product is zero -> divide-by-zero -> NaN), so every atom collapses to
the origin and the ARCSpecies build for atom mapping raises a colliding-atoms
SpeciesError, which aborted the combine and yielded no TS guesses.

That species is used only for connectivity-based atom mapping, so on
SpeciesError we rebuild it from a copy of the zmat with collinear angles
nudged off exactly 0/180 deg (connectivity preserved). The returned TS
coordinates are unaffected: they are built from the original, un-perturbed
zmat_2. Non-linear behavior is byte-for-byte unchanged (the fallback only
runs when the direct build fails).

Adds a regression test exercising a cumulene (O=C=C) acceptor.
…n seed build fails

When CREST's heuristic Z-matrix seed construction yields no seed guesses
(e.g. a linear/cumulene reactive center such as HCCO in H_Abstraction that the
heuristic builder cannot assemble), CREST previously contributed nothing even
when another adapter (e.g. AutoTST) had already produced valid TS guesses for
the same TS.

Add a graceful, opt-in fallback: seed_hub.get_backup_ts_seeds() harvests
successful, non-CREST TS guesses already present on rxn.ts_species.ts_guesses
(the incore TS adapters run sequentially with CREST last, so earlier guesses are
available in-process) and returns them as CREST seeds with empty metadata. The
existing get_wrapper_constraints() machinery then re-derives the family
reactive-atom / distance constraints from the seed geometry itself, so an
externally-built guess is a valid seed without the failed geometry construction.

CrestAdapter.execute_incore() only takes this path when its own seed build
returns nothing; the normal (seed-succeeds) path is unchanged. A feedback-loop
guard (exclude_method='crest') ensures CREST is never seeded from a prior CREST
result, and the pre-existing tiny-system / family gates still apply upstream.

Adds two focused unit tests to crest_test.py covering seed selection
(feedback/failed-guess exclusion, opt_xyz preference, dedup) and the
seed-injection -> constraint-derivation -> input-generation path.
if param.split('_')[0] in ('A', 'AX'):
if value <= delta:
new_vars[param] = delta
elif value >= 180.0 - delta:

Check warning

Code scanning / CodeQL

Redundant comparison Warning

Test is always true, because of
this condition
.

Copilot Autofix

AI 2 days ago

General fix: remove ambiguity in chained range comparisons by computing explicit bounds from validated inputs, then compare against those bounds. This avoids analyzer confusion and makes intent clear.

Best fix here: in arc/job/adapters/ts/heuristics.py, inside _perturb_collinear_zmat_angles, normalize delta to a safe perturbation magnitude in (0, 90) and compute lower_bound/upper_bound once. Then replace direct comparisons against delta / 180.0 - delta with comparisons against these bounds. This keeps current functionality for normal calls (delta=5.0) and prevents pathological values from creating overlapping conditions.

Changes needed:

  • Edit only _perturb_collinear_zmat_angles body (around lines 729–735).
  • No new imports, methods, or external dependencies.
Suggested changeset 1
arc/job/adapters/ts/heuristics.py

Autofix patch

Autofix patch
Run the following command in your local git repository to apply this patch
cat << 'EOF' | git apply
diff --git a/arc/job/adapters/ts/heuristics.py b/arc/job/adapters/ts/heuristics.py
--- a/arc/job/adapters/ts/heuristics.py
+++ b/arc/job/adapters/ts/heuristics.py
@@ -726,13 +726,19 @@
     Returns:
         dict: A copy of ``zmat`` with collinear bond angles perturbed.
     """
+    perturb = abs(float(delta))
+    if perturb >= 90.0:
+        perturb = 89.999
+    lower_bound = perturb
+    upper_bound = 180.0 - perturb
+
     new_vars = dict(zmat['vars'])
     for param, value in zmat['vars'].items():
         if param.split('_')[0] in ('A', 'AX'):
-            if value <= delta:
-                new_vars[param] = delta
-            elif value >= 180.0 - delta:
-                new_vars[param] = 180.0 - delta
+            if value <= lower_bound:
+                new_vars[param] = lower_bound
+            elif value >= upper_bound:
+                new_vars[param] = upper_bound
     return {**zmat, 'vars': new_vars}
 
 
EOF
@@ -726,13 +726,19 @@
Returns:
dict: A copy of ``zmat`` with collinear bond angles perturbed.
"""
perturb = abs(float(delta))
if perturb >= 90.0:
perturb = 89.999
lower_bound = perturb
upper_bound = 180.0 - perturb

new_vars = dict(zmat['vars'])
for param, value in zmat['vars'].items():
if param.split('_')[0] in ('A', 'AX'):
if value <= delta:
new_vars[param] = delta
elif value >= 180.0 - delta:
new_vars[param] = 180.0 - delta
if value <= lower_bound:
new_vars[param] = lower_bound
elif value >= upper_bound:
new_vars[param] = upper_bound
return {**zmat, 'vars': new_vars}


Copilot is powered by AI and may make mistakes. Always verify output.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants