From 7761c2f651d6cf76ce0a490d5da8219ec1beba96 Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Mon, 22 Jun 2026 10:08:39 +0200
Subject: [PATCH 1/9] chore: open PR for issue 229
From b08c5edc443aa14b3723808b6660c4463e287565 Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Mon, 22 Jun 2026 11:17:00 +0200
Subject: [PATCH 2/9] Add paired-state MECP/MECI penalty scanner
Add PairedStateGradientScanner evaluating two electronic surfaces per geometry
from one SCF + one ADC (or LazyMp + ADC for MECP), with joint distinct root
selection and a distinctness guard keeping the tracked roots apart near the
seam. Add a pure Levine-Coe-Martinez penalty objective (mecp_penalty) and an
MECPObjective that drives geomeTRIC through the existing as_pyscf_method
bridge without derivative couplings. Include docs, an example, and unit plus
end-to-end tests; the penalty is cross-checked against geomeTRIC's own
ConicalIntersection engine.
---
adcc/__init__.py | 7 +-
adcc/gradients/__init__.py | 6 +
adcc/gradients/mecp.py | 205 ++++++
adcc/gradients/paired_scanner.py | 590 +++++++++++++++++
adcc/tests/functionality_geomopt_mecp_test.py | 171 +++++
adcc/tests/geomopt_mecp_test.py | 614 ++++++++++++++++++
docs/gradients.rst | 61 +-
examples/optimization/pyscf_adcc_mecp.py | 101 +++
8 files changed, 1751 insertions(+), 4 deletions(-)
create mode 100644 adcc/gradients/mecp.py
create mode 100644 adcc/gradients/paired_scanner.py
create mode 100644 adcc/tests/functionality_geomopt_mecp_test.py
create mode 100644 adcc/tests/geomopt_mecp_test.py
create mode 100644 examples/optimization/pyscf_adcc_mecp.py
diff --git a/adcc/__init__.py b/adcc/__init__.py
index c91d3139..3ad5ae17 100644
--- a/adcc/__init__.py
+++ b/adcc/__init__.py
@@ -47,7 +47,10 @@
from .TwoParticleOperator import TwoParticleOperator
from .TwoParticleDensity import TwoParticleDensity
from .opt_einsum_integration import register_with_opt_einsum
-from .gradients import NuclearGradientScanner, nuclear_gradient
+from .gradients import (NuclearGradientScanner, nuclear_gradient,
+ PairedStateGradientScanner,
+ PairedExcitedStateTarget, PairedGroundExcitedStateTarget,
+ mecp_penalty, MECPObjective)
# This has to be the last set of import
from .guess import (guess_symmetries, guess_zero, guesses_any, guesses_singlet,
@@ -70,6 +73,8 @@
"adc0", "cis", "adc1", "adc2", "adc2x", "adc3",
"cvs_adc0", "cvs_adc1", "cvs_adc2", "cvs_adc2x", "cvs_adc3",
"nuclear_gradient", "NuclearGradientScanner",
+ "PairedStateGradientScanner", "PairedExcitedStateTarget",
+ "PairedGroundExcitedStateTarget", "mecp_penalty", "MECPObjective",
"banner"]
__version__ = "0.18.0"
diff --git a/adcc/gradients/__init__.py b/adcc/gradients/__init__.py
index 8465a150..870ba5f6 100644
--- a/adcc/gradients/__init__.py
+++ b/adcc/gradients/__init__.py
@@ -38,10 +38,16 @@
from .amplitude_response import amplitude_relaxed_densities
from .scanner import (ExcitedStateTarget, GroundStateTarget,
NuclearGradientScanner, density_overlap_score)
+from .paired_scanner import (PairedExcitedStateTarget,
+ PairedGroundExcitedStateTarget,
+ PairedStateGradientScanner)
+from .mecp import mecp_penalty, MECPObjective
__all__ = [
"nuclear_gradient", "NuclearGradientScanner",
"GroundStateTarget", "ExcitedStateTarget", "density_overlap_score",
+ "PairedStateGradientScanner", "PairedExcitedStateTarget",
+ "PairedGroundExcitedStateTarget", "mecp_penalty", "MECPObjective",
]
diff --git a/adcc/gradients/mecp.py b/adcc/gradients/mecp.py
new file mode 100644
index 00000000..bfdb7659
--- /dev/null
+++ b/adcc/gradients/mecp.py
@@ -0,0 +1,205 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2026 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see .
+##
+## ---------------------------------------------------------------------
+"""Penalty-function objective for MECP/MECI optimisations.
+
+The penalty formulation implemented here is the smoothed Levine--Coe--Martinez
+penalty (Levine, Coe and Martinez, *J. Phys. Chem. B* 112, 405, 2008), the same
+form used by geomeTRIC's built-in :class:`ConicalIntersection` engine. It only
+needs the two state energies and gradients -- **no derivative couplings** are
+required, which matches what a :class:`PairedStateGradientScanner` supplies per
+geometry.
+
+This module imports only numpy so that the penalty math stays unit-testable with
+controlled energies/gradients and free of the optional geomeTRIC dependency.
+: class:`MECPObjective` wraps a paired scanner and a penalty into a single
+``(energy, gradient)`` callable plus a ``calc_new`` honouring geomeTRIC's
+flattened custom-engine dict contract; the objective is then driven through
+PySCF's ``as_pyscf_method`` / ``geometric_solver.optimize`` bridge exactly like
+the single-surface scanner.
+"""
+
+from __future__ import annotations
+
+from typing import Callable, Optional
+
+import numpy as np
+
+__all__ = ["mecp_penalty", "MECPObjective"]
+
+# Number of coupled surfaces (MECI/MECP is always a pair).
+_N_STATES = 2
+# geomeTRIC's ConicalIntersection counts each unique state pair once and
+# normalises by ``n_states2 = n_states * (n_states - 1) / 2``; for the two-state
+# MECP/MECI case that is 1. Keep it identical so the oracle cross-check passes.
+_N_STATES2 = _N_STATES * (_N_STATES - 1) // 2
+
+# Default Levine--Coe--Martinez penalty parameters (geomeTRIC's defaults for
+# ``--meci_sigma`` / ``--meci_alpha``). ``alpha`` smooths the otherwise
+# singular penalty at the crossing seam.
+DEFAULT_SIGMA = 3.5
+DEFAULT_ALPHA = 0.025
+
+
+def mecp_penalty(e_lower, g_lower, e_upper, g_upper, *, sigma=DEFAULT_SIGMA,
+ alpha=DEFAULT_ALPHA):
+ """Levine--Coe--Martinez smoothed penalty objective for two surfaces.
+
+ Combines two electronic surfaces ``(e_lower, g_lower)`` and
+ ``(e_upper, g_upper)`` (energies in Hartree, gradients in Hartree/Bohr, any
+ shared shape) into a single penalty objective ``(E, G)`` of the same units
+ suitable for a gradient-based minisation of the crossing seam::
+
+ E_dif = e_upper - e_lower (>= 0 by construction)
+ E_avg = (e_upper + e_lower) / 2
+ E_pen = sigma * E_dif^2 / ((E_dif + alpha) * n_states2)
+ E_obj = E_avg + E_pen
+
+ G_dif = g_upper - g_lower
+ G_avg = (g_upper + g_lower) / 2
+ G_pen = sigma * (E_dif^2 + 2*alpha*E_dif) /
+ ((E_dif + alpha)^2 * n_states2) * G_dif
+ G_obj = G_avg + G_pen
+
+ with ``n_states2 = 1`` for the two-state MECP/MECI case (geomeTRIC counts
+ each unique pair once). The ``alpha``
+ parameter smooths the penalty so the objective stays differentiable at the
+ seam (``E_dif = 0``); the ``sigma`` weight controls how hard the degeneracy
+ is enforced versus minimising the average energy.
+
+ At exact degeneracy (``e_upper == e_lower``) the penalty and its gradient
+ vanish, leaving ``E_obj = E_avg`` and ``G_obj = G_avg`` -- the objective
+ reduces to an average-surface optimisation pinned to the crossing.
+
+ Parameters
+ ----------
+ e_lower, e_upper : float
+ The two state energies; ``e_lower <= e_upper`` is *not* required (the
+ energy difference is formed as ``e_upper - e_lower`` and only its
+ squared phase enters the penalty).
+ g_lower, g_upper : numpy.ndarray
+ The matching state gradients, broadcastable to each other.
+ sigma, alpha : float
+ Penalty strength and smoothing parameter.
+
+ Returns
+ -------
+ (E_obj, G_obj) : (float, numpy.ndarray)
+ The scalar penalty objective and its gradient, with ``G_obj`` shaped
+ like the inputs.
+ """
+ if alpha < 0.0:
+ raise ValueError(f"alpha must be non-negative, got {alpha}.")
+ if sigma < 0.0:
+ raise ValueError(f"sigma must be non-negative, got {sigma}.")
+
+ e_lower = float(e_lower)
+ e_upper = float(e_upper)
+ g_lower = np.asarray(g_lower, dtype=float)
+ g_upper = np.asarray(g_upper, dtype=float)
+
+ e_dif = e_upper - e_lower
+ g_dif = g_upper - g_lower
+ e_avg = 0.5 * (e_upper + e_lower)
+ g_avg = 0.5 * (g_lower + g_upper)
+
+ # Smoothed penalty (Levine--Coe--Martinez / geomeTRIC ConicalIntersection).
+ # At e_dif == 0 the penalty term and its gradient vanish by construction.
+ denom = e_dif + alpha
+ e_pen = sigma * (e_dif * e_dif) / (denom * _N_STATES2)
+ g_pen_scale = (
+ sigma * (e_dif * e_dif + 2.0 * alpha * e_dif)
+ / (denom * denom * _N_STATES2)
+ )
+ g_pen = g_pen_scale * g_dif
+
+ e_obj = e_avg + e_pen
+ g_obj = g_avg + g_pen
+ return float(e_obj), np.asarray(g_obj, dtype=float)
+
+
+class MECPObjective:
+ """Penalty-objective wrapper driving geomeTRIC from a paired scanner.
+
+ Wrap a :class:`PairedStateGradientScanner` (which returns two
+ ``(energy, gradient)`` pairs per geometry) together with the
+ :func:`mecp_penalty` into a single ``__call__`` returning one
+ ``(energy, gradient)`` for PySCF's ``as_pyscf_method`` / geomeTRIC bridge,
+ and a ``calc_new`` returning geomeTRIC's ``{"energy", "gradient"}`` dict
+ with a *flattened* gradient in Hartree/Bohr.
+
+ Parameters
+ ----------
+ scanner : PairedStateGradientScanner
+ The paired scanner supplying both surfaces per geometry.
+ sigma, alpha : float
+ Penalty parameters forwarded to :func:`mecp_penalty`.
+ penalty : callable, optional
+ Custom two-surface penalty with the same signature as
+ :func:`mecp_penalty`; defaults to the Levine--Coe--Martinez smoothed
+ penalty. Useful for injecting a raw energy-difference mode (set
+ ``alpha=0``) or for testing alternative formulations.
+
+ Examples
+ --------
+ >>> objective = MECPObjective(scanner)
+ >>> def energy_and_gradient(mol):
+ ... return objective(mol.atom_coords(unit="Bohr"))
+ >>> method = as_pyscf_method(mol, energy_and_gradient)
+ >>> mol_ci = geometric_solver.optimize(method, maxsteps=20)
+ """
+
+ def __init__(self, scanner, *, sigma=DEFAULT_SIGMA, alpha=DEFAULT_ALPHA,
+ penalty: Optional[Callable] = None):
+ self.scanner = scanner
+ self.sigma = sigma
+ self.alpha = alpha
+ self.penalty = penalty if penalty is not None else mecp_penalty
+ # Bookkeeping of the most recent evaluation for inspection / tests.
+ self.last_energy: Optional[float] = None
+ self.last_gradient = None
+ self.last_pair = None # ((e_lower, g_lower), (e_upper, g_upper))
+
+ def __call__(self, coords):
+ """Return ``(energy, gradient)`` for Cartesian coordinates in Bohr."""
+ (e_lower, g_lower), (e_upper, g_upper) = self.scanner(coords)
+ self.last_pair = ((e_lower, g_lower), (e_upper, g_upper))
+ e_obj, g_obj = self.penalty(
+ e_lower, g_lower, e_upper, g_upper,
+ sigma=self.sigma, alpha=self.alpha,
+ )
+ self.last_energy = e_obj
+ self.last_gradient = np.asarray(g_obj)
+ return e_obj, np.asarray(g_obj)
+
+ def calc_new(self, coords):
+ """geomeTRIC custom-engine entry point.
+
+ ``coords`` is a flattened Cartesian coordinate array in Bohr; the
+ returned gradient is flattened in Hartree/Bohr, matching the contract
+ expected by geomeTRIC's internal engine.
+ """
+ e_obj, g_obj = self(coords)
+ return {
+ "energy": float(e_obj),
+ "gradient": np.asarray(g_obj).ravel(),
+ }
diff --git a/adcc/gradients/paired_scanner.py b/adcc/gradients/paired_scanner.py
new file mode 100644
index 00000000..611ee42a
--- /dev/null
+++ b/adcc/gradients/paired_scanner.py
@@ -0,0 +1,590 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2026 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see .
+##
+## ---------------------------------------------------------------------
+"""Paired-surface nuclear-gradient scanner for MECP/MECI optimisations.
+
+: class:`PairedStateGradientScanner` evaluates *two* electronic surfaces at a
+single geometry from one SCF + one ADC (or, for MECP, one ``LazyMp`` plus one
+ADC) and returns both ``(energy, gradient)`` pairs. It reuses the
+single-surface : class:`NuclearGradientScanner` plumbing (SCF lifecycle, AO
+density-overlap root tracking, coordinate handling) and adds a **joint distinct
+root selection** with a distinctness guard so the two tracked roots never
+collapse onto the same candidate state.
+
+A pure penalty objective that combines the two surfaces into one
+``(energy, gradient)`` for a geomeTRIC-driven optimisation lives in a separate
+module; this file owns only the paired scanner and its target dataclasses.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+from typing import Any, Optional, Sequence
+
+import numpy as np
+
+from adcc.LazyMp import LazyMp
+from adcc.Excitation import Excitation
+
+from adcc.gradients.scanner import (
+ GroundStateTarget,
+ NuclearGradientScanner,
+ TrackingResult,
+ _TrackingDescriptor,
+ _check_ground_state_level,
+ _mp_level,
+ density_overlap_score,
+)
+
+
+__all__ = [
+ "PairedExcitedStateTarget",
+ "PairedGroundExcitedStateTarget",
+ "PairedStateGradientScanner",
+ "density_overlap_score",
+]
+
+
+# ---------------------------------------------------------------------------
+# Paired target dataclasses
+# ---------------------------------------------------------------------------
+
+@dataclass(frozen=True)
+class PairedExcitedStateTarget:
+ """Two excited-state ADC roots evaluated from a single ``run_adc`` call.
+
+ Used for MECI scans where both surfaces come from the same excited-state
+ manifold. ``state_indices`` are the *positional* seeds for the two tracked
+ roots; overlap tracking (when enabled) follows their state character across
+ geometries while keeping the two roots distinct.
+ """
+
+ method: str
+ state_indices: tuple[int, int]
+ run_adc_kwargs: dict[str, Any] = field(default_factory=dict)
+
+ def __post_init__(self):
+ if len(self.state_indices) != 2:
+ raise ValueError(
+ "PairedExcitedStateTarget.state_indices must contain exactly two "
+ f"indices, got {self.state_indices!r}."
+ )
+ # Immutable tuple of plain ints.
+ object.__setattr__(
+ self, "state_indices",
+ (int(self.state_indices[0]), int(self.state_indices[1])),
+ )
+
+ def kwargs(self) -> dict[str, Any]:
+ """Return keyword arguments forwarded to : func:`adcc.run_adc`."""
+ ret = dict(self.run_adc_kwargs)
+ ret["method"] = self.method
+ return ret
+
+
+@dataclass(frozen=True)
+class PairedGroundExcitedStateTarget:
+ """One ground-state MP surface plus one excited-state ADC surface (MECP).
+
+ The ground surface is always MP2 (the only level for which adcc provides a
+ nuclear gradient) and is evaluated via : class:`adcc.LazyMp` /
+ ``LazyMp.energy(2)``. The excited surface is the single ADC root at
+ ``state_index`` of the same SCF reference. Only the excited root is
+ overlap-tracked; the ground side is uniquely defined by the SCF reference.
+ """
+
+ method: str
+ state_index: int = 0
+ run_adc_kwargs: dict[str, Any] = field(default_factory=dict)
+ level: int = 2
+
+ def __post_init__(self):
+ _check_ground_state_level(self.level)
+ object.__setattr__(self, "state_index", int(self.state_index))
+
+ def kwargs(self) -> dict[str, Any]:
+ """Return keyword arguments forwarded to : func:`adcc.run_adc`."""
+ ret = dict(self.run_adc_kwargs)
+ ret["method"] = self.method
+ return ret
+
+
+# ---------------------------------------------------------------------------
+# Paired scanner
+# ---------------------------------------------------------------------------
+
+class PairedStateGradientScanner(NuclearGradientScanner):
+ """Callable scanner returning two ``(energy, gradient)`` surface pairs.
+
+ ``scfres`` is a configured PySCF SCF object; its molecule and SCF settings
+ are the template for every scanner call and only the geometry changes.
+ Coordinates are Cartesian in Bohr, energies are Hartree and gradients are
+ Hartree/Bohr.
+
+ Two normalised target forms are supported:
+
+ * **MECI** -- two excited roots from one ADC solve, via
+ ``PairedExcitedStateTarget`` or the ``method=..., states=(i, j)``
+ convenience.
+ * **MECP** -- one MP2 ground surface plus one excited ADC root, via a
+ ``PairedGroundExcitedStateTarget`` or the ``lower="mp2", upper=k``
+ convenience (``method`` supplies the excited-state ADC method).
+
+ ``__call__`` returns ``((e_lower, g_lower), (e_upper, g_upper))`` with the
+ two surfaces sorted by energy so that ``e_lower <= e_upper``. Both surfaces
+ are also kept in *slot* order on :attr:`last_energies` / :attr:`last_gradients`
+ (slot 0 = first tracked root / ground, slot 1 = second tracked root /
+ excited) for diagnostics.
+ """
+
+ def __init__(self, scfres, *,
+ target: Optional[Any] = None,
+ method: Optional[str] = None,
+ states: Optional[Sequence[int]] = None,
+ state_indices: Optional[Sequence[int]] = None,
+ lower: Optional[Any] = None,
+ upper: Optional[int] = None,
+ mp_level: int = 2,
+ follow: str = "overlap",
+ tracking_min_score: float = 0.0,
+ tracking_min_gap: float = 0.0,
+ gradient_kwargs: Optional[dict[str, Any]] = None,
+ **run_adc_kwargs):
+ # Reuse the single-surface plumbing: PySCF import/validation, the SCF
+ # scanner lifecycle, coordinate helpers and the gradient kwargs store.
+ # The parent normalises a placeholder single-surface target (an MP2
+ # ground-state target by default); the paired target is built below and
+ # takes precedence over the placeholder on every code path.
+ super().__init__(
+ scfres,
+ follow=follow,
+ tracking_min_score=tracking_min_score,
+ tracking_min_gap=tracking_min_gap,
+ gradient_kwargs=gradient_kwargs,
+ mp_level=2,
+ )
+
+ self.paired_target = self._normalise_paired_target(
+ target, method, states, state_indices, lower, upper, mp_level,
+ run_adc_kwargs,
+ )
+
+ # Plural tracking state. For MECP slot 0 (ground) is never tracked:
+ # its descriptor / index stay ``None`` forever.
+ self.previous_descriptors: tuple[Optional[_TrackingDescriptor],
+ Optional[_TrackingDescriptor]] = (
+ None, None)
+ self.previous_indices: tuple[Optional[int], Optional[int]] = (None, None)
+
+ # Last-call diagnostics (slot order).
+ self.last_scf = None
+ self.previous_scf = None
+ self.last_states = None
+ self.last_targets: tuple[Optional[Any], Optional[Any]] = (None, None)
+ self.last_excitations: tuple[Optional[Excitation],
+ Optional[Excitation]] = (None, None)
+ self.last_gradients: tuple[np.ndarray, np.ndarray] = (None, None)
+ self.last_gradient_results: tuple[Any, Any] = (None, None)
+ self.last_energies: tuple[float, float] = (None, None)
+ self.last_trackings: tuple[Optional[TrackingResult],
+ Optional[TrackingResult]] = (None, None)
+ # Energy-sorted insertion order of the two surfaces in the last call:
+ # ``last_pair_order = (lo_slot, hi_slot)``.
+ self.last_pair_order: tuple[int, int] = (0, 1)
+
+ # -- public API --------------------------------------------------------
+
+ def __call__(self, coords):
+ """Return ``((e_lower, g_lower), (e_upper, g_upper))`` for ``coords``.
+
+ ``coords`` are Cartesian coordinates in Bohr with shape
+ ``(natoms, 3)`` or ``(3 * natoms,)``. The two surfaces are sorted by
+ energy so the returned lower energy is never larger than the upper one.
+ """
+ coords = self._coords_array(coords)
+ scfres = self._run_scf(coords)
+ targets = self._build_target(scfres)
+
+ from adcc.gradients import nuclear_gradient
+ grad_results = [
+ nuclear_gradient(t, **self.gradient_kwargs) for t in targets
+ ]
+ grads = [np.asarray(g.total) for g in grad_results]
+ energies = []
+ for t in targets:
+ if isinstance(t, LazyMp):
+ energies.append(float(t.energy(self.paired_target.level)))
+ else:
+ energies.append(float(t.total_energy))
+
+ self.last_scf = scfres
+ self.previous_scf = scfres
+ self.last_targets = tuple(targets)
+ self.last_gradient_results = tuple(grad_results)
+ self.last_gradients = tuple(grads)
+ self.last_energies = tuple(energies)
+
+ order = np.argsort(energies)
+ lo, hi = int(order[0]), int(order[1])
+ self.last_pair_order = (lo, hi)
+ return ((energies[lo], grads[lo]), (energies[hi], grads[hi]))
+
+ def calc_new(self, coords):
+ """geomeTRIC-style entry point returning *both* surface pairs.
+
+ Unlike the single-surface scanner, this returns a dictionary carrying
+ the two energies and (flattened) gradients separately rather than a
+ single ``energy``/``gradient`` pair. A pure penalty objective combining
+ the two surfaces into one geomeTRIC-bound surface lives in a separate
+ module and drives the optimizer; this method is intended for inspection
+ and testing.
+ """
+ (e_lo, g_lo), (e_hi, g_hi) = self(coords)
+ g_lo = np.asarray(g_lo)
+ g_hi = np.asarray(g_hi)
+ return {
+ "energies": np.array([e_lo, e_hi]),
+ "gradients": np.stack([g_lo, g_hi]),
+ "energy_lower": float(e_lo),
+ "energy_upper": float(e_hi),
+ "gradient_lower": g_lo.ravel(),
+ "gradient_upper": g_hi.ravel(),
+ }
+
+ # -- target normalisation ---------------------------------------------
+
+ def _normalise_paired_target(self, target, method, states, state_indices,
+ lower, upper, mp_level, run_adc_kwargs):
+ if isinstance(target, PairedExcitedStateTarget):
+ return target
+ if isinstance(target, PairedGroundExcitedStateTarget):
+ # ``__post_init__`` already validated the MP level.
+ return target
+ if isinstance(target, str):
+ method = target
+ run_adc_kwargs = dict(run_adc_kwargs)
+
+ states = states if states is not None else state_indices
+
+ # MECP form: lower=, upper=.
+ if lower is not None and upper is not None:
+ if isinstance(lower, GroundStateTarget):
+ level = lower.level
+ elif isinstance(lower, str):
+ level = _mp_level(lower, mp_level)
+ elif isinstance(lower, int):
+ level = lower
+ else:
+ raise TypeError(
+ "lower must be 'mp2', an int MP level, or a "
+ "GroundStateTarget, got "
+ f"{type(lower).__name__}."
+ )
+ _check_ground_state_level(level)
+ if not isinstance(upper, int):
+ raise TypeError(
+ "upper must be an int excited-state index, got "
+ f"{type(upper).__name__}."
+ )
+ if method is None:
+ raise ValueError(
+ "An ADC method is required for the excited side of an "
+ "MECP scanner (pass method=...)."
+ )
+ return PairedGroundExcitedStateTarget(
+ method=method, state_index=upper,
+ run_adc_kwargs=dict(run_adc_kwargs), level=level,
+ )
+
+ # MECI form: two excited indices from one ADC solve.
+ if states is not None:
+ if not isinstance(states, (tuple, list)):
+ raise TypeError(
+ "states must be a 2-tuple of excited-state indices, got "
+ f"{type(states).__name__}."
+ )
+ states = tuple(states)
+ if len(states) != 2:
+ raise ValueError(
+ "states must contain exactly two excited-state indices, "
+ f"got {len(states)}."
+ )
+ if method is None:
+ raise ValueError(
+ "An ADC method is required for a paired excited-state "
+ "scanner (pass method=...)."
+ )
+ return PairedExcitedStateTarget(
+ method=method, state_indices=tuple(states),
+ run_adc_kwargs=dict(run_adc_kwargs),
+ )
+
+ raise ValueError(
+ "PairedStateGradientScanner needs a paired target: pass a "
+ "PairedExcitedStateTarget / PairedGroundExcitedStateTarget, or the "
+ "convenience forms method=..., states=(i, j) (MECI) / "
+ "lower='mp2', upper=k (MECP)."
+ )
+
+ # -- per-geometry build -----------------------------------------------
+
+ def _build_target(self, scfres):
+ """Run SCF-anchored adcc once and return the two surface targets."""
+ from adcc import ReferenceState, run_adc
+
+ if isinstance(self.paired_target, PairedExcitedStateTarget):
+ states = run_adc(scfres, **self.paired_target.kwargs())
+ self.last_states = states
+ chosen, descriptors, trackings = self._select_excitations(
+ states, scfres.mol
+ )
+ self.last_trackings = tuple(trackings)
+ self.last_excitations = (
+ states.excitations[chosen[0]], states.excitations[chosen[1]],
+ )
+ self.last_targets = self.last_excitations
+ self.previous_descriptors = (
+ descriptors[chosen[0]], descriptors[chosen[1]],
+ )
+ self.previous_indices = (chosen[0], chosen[1])
+ return self.last_targets
+
+ # MECP: MP2 ground + single tracked excited root from one SCF + one ADC.
+ mp = LazyMp(ReferenceState(scfres))
+ states = run_adc(scfres, **self.paired_target.kwargs())
+ self.last_states = states
+ excitation, tracking, descriptor, idx = self._select_excited_root(
+ states, scfres.mol,
+ self.previous_descriptors[1], self.previous_indices[1],
+ self.paired_target.state_index,
+ )
+ self.last_trackings = (None, tracking)
+ self.last_excitations = (None, excitation)
+ self.last_targets = (mp, excitation)
+ self.previous_descriptors = (None, descriptor)
+ self.previous_indices = (None, idx)
+ return self.last_targets
+
+ # -- root selection ----------------------------------------------------
+
+ def _select_excitations(self, states, mol):
+ """Jointly select two distinct excited roots for the MECI pair.
+
+ Reuses the single-surface AO density-overlap tracking machinery via
+ :meth:`NuclearGradientScanner._tracking_score` and adds a joint
+ distinctness guard so the two tracked roots never collapse onto the same
+ candidate state. Returns ``(chosen, descriptors, trackings)`` where
+ ``chosen`` is the two selected positional indices (slot order).
+ """
+ excitations = states.excitations
+ if not excitations:
+ raise RuntimeError(
+ "ADC calculation did not return any excited states."
+ )
+ if len(excitations) < 2:
+ raise RuntimeError(
+ "Paired excited-state scanner needs at least two excited "
+ f"states, got {len(excitations)}."
+ )
+ n = len(excitations)
+ descriptors = [
+ self._descriptor(excitation, mol) for excitation in excitations
+ ]
+
+ unavailable: set[str] = set()
+ score_vectors: list[Optional[np.ndarray]] = [None, None]
+ for slot in range(2):
+ prev = self.previous_descriptors[slot]
+ if self.follow == "overlap" and prev is not None:
+ score_vectors[slot] = np.array([
+ self._tracking_score(prev, descriptors[k], mol, unavailable)
+ for k in range(n)
+ ])
+
+ # Slots without a score vector fall back to their positional seed index.
+ target_indices = tuple(self.paired_target.state_indices)
+ fixed = [None, None]
+ for slot in range(2):
+ if score_vectors[slot] is None:
+ idx = int(target_indices[slot])
+ if idx >= n or idx < -n:
+ raise ValueError(
+ f"state_index {idx} is out of range for {n} computed "
+ "states."
+ )
+ fixed[slot] = idx % n
+
+ free_slots = [s for s in range(2) if score_vectors[s] is not None]
+
+ if len(free_slots) == 0:
+ chosen = [fixed[0], fixed[1]]
+ elif len(free_slots) == 1:
+ s = free_slots[0]
+ other = 1 - s
+ other_idx = fixed[other]
+ order = np.argsort(score_vectors[s])[::-1]
+ pick = None
+ for o in order:
+ if int(o) != other_idx:
+ pick = int(o)
+ break
+ if pick is None:
+ raise RuntimeError(
+ "Root tracking could not find two distinct excited roots "
+ "for the paired scanner."
+ )
+ chosen = [None, None]
+ chosen[s] = pick
+ chosen[other] = other_idx
+ else:
+ # Both slots track: jointly maximise the summed overlap score over
+ # all distinct ordered pairs. Near-degenerate ties (combined scores
+ # within ``eps``) are broken towards the pair with the smaller
+ # excitation-energy gap, i.e. the seam the optimisation is hunting.
+ used = {f for f in fixed if f is not None}
+ available = [k for k in range(n) if k not in used]
+ if len(available) < 2:
+ raise RuntimeError(
+ "Not enough distinct excited states to keep the two "
+ "tracked roots apart."
+ )
+ eps = 1e-9
+ best_pair = None
+ best_score = -np.inf
+ best_tie_gap = np.inf
+ for i in available:
+ for j in available:
+ if i == j:
+ continue
+ total = float(score_vectors[0][i] + score_vectors[1][j])
+ tie_gap = abs(
+ float(excitations[i].excitation_energy)
+ - float(excitations[j].excitation_energy)
+ )
+ better = total > best_score + eps
+ tie = abs(total - best_score) <= eps and tie_gap < best_tie_gap
+ if better or tie:
+ best_pair = (i, j)
+ best_score = total
+ best_tie_gap = tie_gap
+ chosen = [best_pair[0], best_pair[1]]
+
+ if chosen[0] == chosen[1]:
+ raise RuntimeError(
+ "Distinctness guard: the two tracked roots collapsed onto "
+ f"state {chosen[0]}. Cannot keep the paired surfaces apart."
+ )
+
+ trackings: list[Optional[TrackingResult]] = [None, None]
+ for slot in range(2):
+ sv = score_vectors[slot]
+ idx = chosen[slot]
+ if sv is None:
+ # Positional seeding: mirror the single-surface scanner, which
+ # records no tracking diagnostic on the seeding call.
+ continue
+ other = chosen[1 - slot]
+ # Gap to the best *alternative* root available to this slot, i.e.
+ # excluding the partner's selected root (the distinctness guard
+ # means there is always at least one other candidate here).
+ remaining = [k for k in range(n) if k != idx and k != other]
+ if remaining:
+ second = max(remaining, key=lambda k: sv[k])
+ gap = float(sv[idx] - sv[second])
+ else:
+ gap = float("inf")
+ prev_idx = self.previous_indices[slot]
+ switched = prev_idx is not None and idx != prev_idx
+ trackings[slot] = TrackingResult(
+ index=idx, scores=sv, best_score=float(sv[idx]), gap=gap,
+ previous_index=prev_idx, switched=bool(switched),
+ unavailable_channels=tuple(sorted(unavailable)),
+ )
+ if sv[idx] < self.tracking_min_score:
+ raise RuntimeError(
+ f"Root tracking failed for paired slot {slot}: best "
+ f"state-character overlap {sv[idx]:.6g} is below "
+ f"threshold {self.tracking_min_score:.6g}."
+ )
+ if gap < self.tracking_min_gap:
+ raise RuntimeError(
+ f"Root tracking is ambiguous for paired slot {slot}: "
+ f"best and second-best overlaps differ by {gap:.6g}."
+ )
+ return chosen, descriptors, trackings
+
+ def _select_excited_root(self, states, mol, prev_descriptor, prev_index,
+ seed_index):
+ """Select a single tracked excited root (MECP excited side).
+
+ A faithful paired-side copy of
+ :meth:`NuclearGradientScanner._select_excitation` that keeps the
+ previous-descriptor / previous-index state local to this slot. Returns
+ ``(excitation, tracking, descriptor, index)``.
+ """
+ excitations = states.excitations
+ if not excitations:
+ raise RuntimeError(
+ "ADC calculation did not return any excited states."
+ )
+ n = len(excitations)
+
+ if self.follow == "index" or prev_descriptor is None:
+ idx = int(seed_index)
+ if idx >= n or idx < -n:
+ raise ValueError(
+ f"state_index {idx} is out of range for {n} computed states."
+ )
+ idx %= n
+ descriptor = self._descriptor(excitations[idx], mol)
+ return excitations[idx], None, descriptor, idx
+
+ unavailable: set[str] = set()
+ descriptors = [
+ self._descriptor(excitation, mol) for excitation in excitations
+ ]
+ scores = np.array([
+ self._tracking_score(prev_descriptor, descriptors[k], mol, unavailable)
+ for k in range(n)
+ ])
+ order = np.argsort(scores)[::-1]
+ best = int(order[0])
+ gap = (float(scores[best] - scores[int(order[1])])
+ if len(order) > 1 else float("inf"))
+ switched = prev_index is not None and best != prev_index
+ tracking = TrackingResult(
+ index=best, scores=scores, best_score=float(scores[best]), gap=gap,
+ previous_index=prev_index, switched=bool(switched),
+ unavailable_channels=tuple(sorted(unavailable)),
+ )
+ if scores[best] < self.tracking_min_score:
+ raise RuntimeError(
+ "Root tracking failed: best state-character overlap "
+ f"{scores[best]:.6g} is below threshold "
+ f"{self.tracking_min_score:.6g}."
+ )
+ if len(order) > 1 and gap < self.tracking_min_gap:
+ raise RuntimeError(
+ "Root tracking is ambiguous: best and second-best overlaps "
+ f"differ by {gap:.6g}."
+ )
+ return excitations[best], tracking, descriptors[best], best
diff --git a/adcc/tests/functionality_geomopt_mecp_test.py b/adcc/tests/functionality_geomopt_mecp_test.py
new file mode 100644
index 00000000..aa422c46
--- /dev/null
+++ b/adcc/tests/functionality_geomopt_mecp_test.py
@@ -0,0 +1,171 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2026 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see .
+##
+## ---------------------------------------------------------------------
+"""End-to-end MECP/MECI optimisation smoke tests.
+
+These drive the :class:`PairedStateGradientScanner` + :class:`MECPObjective`
+penalty through geomeTRIC (via PySCF's geomopt bridge) and are skipped cleanly
+when either optional dependency is unavailable. They exercise the paired
+root-tracking + distinctness guard across real optimisation steps; the heavy
+full-optimisation convergence is left to the example script and these tests
+keep ``maxsteps`` tiny.
+"""
+import importlib.util
+
+import numpy as np
+import pytest
+
+import adcc
+import adcc.backends
+
+
+def _missing(*modules):
+ return [m for m in modules if importlib.util.find_spec(m) is None]
+
+
+_required = ["pyscf", "geometric"]
+pytestmark = pytest.mark.skipif(
+ "pyscf" not in adcc.backends.available() or _missing(*_required),
+ reason="PySCF and geomeTRIC are required for end-to-end MECP/MECI tests.",
+)
+
+
+def _ethylene_scf():
+ """A small twisted ethylene SCF: the textbook minimal MECI system.
+
+ A near-90 deg torsion about the C=C bond drives the two lowest singlet
+ surfaces toward a conical intersection, keeping the system cheap at
+ STO-3G. This fixture is authored fresh (no ethylene fixture exists yet).
+ """
+ from pyscf import gto, scf
+ # Twist one CH2 group ~80 deg to start close to the crossing seam.
+ tw = np.deg2rad(80.0)
+ c1 = np.array([0.0, 0.0, 0.0])
+ c2 = np.array([1.34, 0.0, 0.0])
+ # Left CH2 (untwisted).
+ h1l = c1 + np.array([0.0, 0.63, 0.0])
+ h2l = c1 + np.array([0.0, -0.63, 0.0])
+ # Right CH2 (twisted about the C=C axis = x).
+
+ def twist(p, angle, pivot):
+ r = p - pivot
+ rot = np.array([[1.0, 0.0, 0.0],
+ [0.0, np.cos(angle), -np.sin(angle)],
+ [0.0, np.sin(angle), np.cos(angle)]])
+ return pivot + rot @ r
+ h1r = twist(c2 + np.array([0.0, 0.63, 0.0]), tw, c2)
+ h2r = twist(c2 + np.array([0.0, -0.63, 0.0]), tw, c2)
+ atoms = ["C", "C", "H", "H", "H", "H"]
+ coords = np.stack([c1, c2, h1l, h2l, h1r, h2r])
+ atom_str = "\n".join(
+ f"{sym} {xyz[0]:.8f} {xyz[1]:.8f} {xyz[2]:.8f}"
+ for sym, xyz in zip(atoms, coords)
+ )
+ mol = gto.M(
+ atom=atom_str, basis="sto-3g", unit="Angstrom",
+ symmetry=False, verbose=0, parse_arg=False,
+ )
+ mf = scf.RHF(mol)
+ mf.conv_tol = 1e-10
+ mf.conv_tol_grad = 1e-7
+ return mf
+
+
+def test_paired_scanner_objective_drives_meci_optimization():
+ from pyscf.geomopt import as_pyscf_method, geometric_solver
+
+ scfres = _ethylene_scf()
+ scanner = adcc.PairedStateGradientScanner(
+ scfres, method="adc2", states=(0, 1), n_singlets=3,
+ follow="overlap", conv_tol=1e-8,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ objective = adcc.MECPObjective(scanner) # default LCM penalty
+
+ seen_gaps = []
+
+ def energy_and_gradient(mol_at_step):
+ energy, gradient = objective(mol_at_step.atom_coords(unit="Bohr"))
+ e_lo, e_hi = objective.last_pair[0][0], objective.last_pair[1][0]
+ seen_gaps.append(abs(e_hi - e_lo))
+ return energy, gradient
+
+ method = as_pyscf_method(scfres.mol, energy_and_gradient)
+ geometric_solver.optimize(
+ method, maxsteps=3,
+ convergence_grms=1e-4, convergence_gmax=2e-4,
+ )
+
+ # Multiple steps ran and both surfaces were evaluated every step.
+ assert len(seen_gaps) >= 2
+ assert all(np.isfinite(gap) for gap in seen_gaps)
+ assert scanner.last_trackings is not None
+ # Tracking across steps kept the two followed roots distinct throughout.
+ slot0, slot1 = scanner.last_trackings
+ assert slot0 is not None and slot1 is not None
+ assert slot0.index != slot1.index
+
+
+def test_mecp_objective_calc_new_contract_end_to_end():
+ # Reuse the ethylene system to confirm the MECPObjective's calc_new returns
+ # the geomeTRIC custom-engine dict with finite energy and a flattened,
+ # Bohr-scaled gradient.
+ scfres = _ethylene_scf()
+ scanner = adcc.PairedStateGradientScanner(
+ scfres, method="adc2", states=(0, 1), n_singlets=3,
+ follow="index", conv_tol=1e-8,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ objective = adcc.MECPObjective(scanner)
+ result = objective.calc_new(scanner.initial_coords.ravel())
+ assert set(result) == {"energy", "gradient"}
+ assert isinstance(result["energy"], float)
+ assert np.isfinite(result["energy"])
+ assert result["gradient"].shape == (3 * scanner.natoms,)
+ assert np.all(np.isfinite(result["gradient"]))
+
+
+def test_mecp_ground_excited_pair_drives_optimization():
+ from pyscf.geomopt import as_pyscf_method, geometric_solver
+
+ scfres = _ethylene_scf()
+ scanner = adcc.PairedStateGradientScanner(
+ scfres, method="adc2", lower="mp2", upper=0, n_singlets=3,
+ follow="index", conv_tol=1e-8,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ objective = adcc.MECPObjective(scanner)
+
+ def energy_and_gradient(mol_at_step):
+ return objective(mol_at_step.atom_coords(unit="Bohr"))
+
+ method = as_pyscf_method(scfres.mol, energy_and_gradient)
+ # A couple of steps confirm the MECP objective drives geomeTRIC without
+ # errors; full seam convergence is not required here.
+ geometric_solver.optimize(
+ method, maxsteps=2,
+ convergence_grms=1e-4, convergence_gmax=2e-4,
+ )
+ # The ground/excited pair stayed finite and ordered.
+ e_lo, e_hi = objective.last_pair[0][0], objective.last_pair[1][0]
+ assert np.isfinite(e_lo) and np.isfinite(e_hi)
+ assert e_lo <= e_hi
diff --git a/adcc/tests/geomopt_mecp_test.py b/adcc/tests/geomopt_mecp_test.py
new file mode 100644
index 00000000..0ee707de
--- /dev/null
+++ b/adcc/tests/geomopt_mecp_test.py
@@ -0,0 +1,614 @@
+#!/usr/bin/env python3
+## vi: tabstop=4 shiftwidth=4 softtabstop=4 expandtab
+## ---------------------------------------------------------------------
+##
+## Copyright (C) 2026 by the adcc authors
+##
+## This file is part of adcc.
+##
+## adcc is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## by the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## adcc is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with adcc. If not, see .
+##
+## ---------------------------------------------------------------------
+"""Fast unit tests for the paired-state MECP/MECI scanner and penalty objective.
+
+These mirror the gating style of ``geomopt_scanner_test.py`` (PySCF-only, no
+geomeTRIC dependency) and cover: the pure penalty math (finite-difference
+gradient consistency, the Levine-Coe-Martinez formula, degeneracy limits and
+parameter validation), the :class:`MECPObjective` geomeTRIC bridge contract, the
+paired-engine mechanics (two distinct surfaces from one SCF + one ADC/MMP), the
+distinctness guard, and paired-target normalisation/validation.
+"""
+import types
+
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+import adcc
+import adcc.backends
+from adcc.gradients.mecp import mecp_penalty, MECPObjective, _N_STATES2
+from adcc.gradients.paired_scanner import (
+ PairedExcitedStateTarget,
+ PairedGroundExcitedStateTarget,
+)
+from adcc.gradients.scanner import _TrackingDescriptor
+
+
+pytestmark = pytest.mark.skipif(
+ "pyscf" not in adcc.backends.available(), reason="PySCF not found."
+)
+
+
+def _h2o_scf():
+ from pyscf import gto, scf
+ mol = gto.M(
+ atom="""
+ O 0 0 0
+ H 0 0 1.795239827225189
+ H 1.693194615993441 0 -0.599043184453037
+ """,
+ basis="sto-3g", unit="Bohr", symmetry=False,
+ verbose=0, parse_arg=False,
+ )
+ mf = scf.RHF(mol)
+ mf.conv_tol = 1e-11
+ mf.conv_tol_grad = 1e-9
+ return mf
+
+
+def _mecp_scf():
+ """A small twisted ethylene-ish scratch SCF with >=2 singlet roots."""
+ from pyscf import gto, scf
+ mol = gto.M(
+ atom="""
+ C 0.0 0.0 0.0
+ C 1.3 0.0 0.0
+ H 0.0 1.0 0.0
+ H 0.0 -1.0 0.0
+ H 1.3 1.0 0.0
+ H 1.3 -1.0 0.0
+ """,
+ basis="sto-3g", unit="Angstrom", symmetry=False,
+ verbose=0, parse_arg=False,
+ )
+ mf = scf.RHF(mol)
+ mf.conv_tol = 1e-11
+ mf.conv_tol_grad = 1e-9
+ return mf
+
+
+# ---------------------------------------------------------------------------
+# Penalty math (controlled energies/gradients) -- no SCF needed.
+# ---------------------------------------------------------------------------
+
+def test_n_states2_is_one_for_two_state_pair():
+ assert _N_STATES2 == 1
+
+
+def test_penalty_gradient_matches_finite_difference():
+ rng = np.random.default_rng(42)
+ e_lo, e_hi = -1.3, -1.1
+ g_lo = rng.normal(size=(2, 3))
+ g_hi = rng.normal(size=(2, 3))
+
+ energy, gradient = mecp_penalty(e_lo, g_lo, e_hi, g_hi,
+ sigma=3.5, alpha=0.025)
+
+ h = 1e-7
+
+ def efun(el, eh):
+ return mecp_penalty(el, g_lo, eh, g_hi, sigma=3.5, alpha=0.025)[0]
+
+ d_lo = (efun(e_lo + h, e_hi) - efun(e_lo - h, e_hi)) / (2 * h)
+ d_hi = (efun(e_lo, e_hi + h) - efun(e_lo, e_hi - h)) / (2 * h)
+ gradient_fd = d_lo * g_lo + d_hi * g_hi
+ assert gradient.shape == g_lo.shape
+ assert_allclose(gradient, gradient_fd, atol=1e-6)
+
+
+def test_penalty_energy_matches_hand_computed_lcm_formula():
+ # EAvg = (-1.3 + -1.1)/2 = -1.2; EDif = 0.2; with sigma=3.5, alpha=0.025,
+ # n_states2=1: EPen = 3.5 * 0.04 / (0.225) = 0.6222..., E = EAvg + EPen.
+ e_lo, e_hi = -1.3, -1.1
+ g_lo = np.zeros((1, 3))
+ g_hi = np.zeros((1, 3))
+ energy, _ = mecp_penalty(e_lo, g_lo, e_hi, g_hi, sigma=3.5, alpha=0.025)
+ e_avg = 0.5 * (e_lo + e_hi)
+ e_dif = e_hi - e_lo
+ e_pen = 3.5 * e_dif ** 2 / ((e_dif + 0.025) * 1)
+ assert energy == pytest.approx(e_avg + e_pen)
+
+
+def test_penalty_penalty_vanishes_at_exact_degeneracy():
+ g_lo = np.array([[1.0, 0.0, 0.0]])
+ g_hi = np.array([[0.0, 2.0, 0.0]])
+ e = -1.2
+ energy, gradient = mecp_penalty(e, g_lo, e, g_hi, sigma=3.5, alpha=0.025)
+ # At EDif == 0 the penalty and its gradient vanish, so the objective reduces
+ # to the average surface.
+ assert energy == pytest.approx(e)
+ assert_allclose(gradient, 0.5 * (g_lo + g_hi))
+
+
+def test_penalty_objective_is_flat_when_both_gradients_zero_at_degeneracy():
+ e = -1.0
+ zero = np.zeros((3, 3))
+ energy, gradient = mecp_penalty(e, zero, e, zero)
+ assert energy == pytest.approx(e)
+ assert_allclose(gradient, zero)
+
+
+def test_penalty_alpha_zero_is_raw_squared_difference_mode():
+ e_lo, e_hi = -1.3, -1.1
+ g_lo = np.zeros((1, 3))
+ g_hi = np.array([[1.0, 0.0, 0.0]])
+ # alpha == 0: the penalty becomes sigma * EDif (energy-difference objective)
+ # so the penalty's energy contribution is sigma * EDif and the gradient
+ # contribution is sigma * GDif about the upper-lower split.
+ energy, gradient = mecp_penalty(e_lo, g_lo, e_hi, g_hi, sigma=2.0, alpha=0.0)
+ e_dif = e_hi - e_lo
+ assert energy == pytest.approx(0.5 * (e_lo + e_hi) + 2.0 * e_dif)
+ assert_allclose(gradient, 0.5 * (g_lo + g_hi) + 2.0 * (g_hi - g_lo))
+
+
+def test_penalty_parameters_are_validated():
+ with pytest.raises(ValueError, match="sigma"):
+ mecp_penalty(0.0, np.zeros(3), 1.0, np.zeros(3), sigma=-1.0)
+ with pytest.raises(ValueError, match="alpha"):
+ mecp_penalty(0.0, np.zeros(3), 1.0, np.zeros(3), alpha=-0.1)
+
+
+def test_penalty_tolerates_unsorted_input_without_sort_inversion():
+ # Passing e_upper < e_lower keeps e_dif = e_upper - e_lower (negative); the
+ # squared terms stay real and finite, the objective stays finite. This only
+ # documents that the function tolerates unsorted input; callers (the paired
+ # scanner) always feed energy-sorted surfaces.
+ energy, gradient = mecp_penalty(-1.1, np.zeros(3), -1.3, np.zeros(3))
+ assert np.isfinite(energy)
+ assert_allclose(gradient, np.zeros(3))
+
+
+# ---------------------------------------------------------------------------
+# MECPObjective bridge contract.
+# ---------------------------------------------------------------------------
+
+def test_mecp_objective_returns_combined_energy_and_gradient():
+ g_lo = np.ones((2, 3))
+ g_hi = 2 * np.ones((2, 3))
+ energy_direct, gradient_direct = mecp_penalty(
+ -1.3, g_lo, -1.1, g_hi, sigma=3.5, alpha=0.025,
+ )
+
+ class _FakeScanner:
+ def __call__(self, coords):
+ return ((-1.3, g_lo), (-1.1, g_hi))
+
+ obj = MECPObjective(_FakeScanner(), sigma=3.5, alpha=0.025)
+ energy, gradient = obj("coords-ignored")
+ assert energy == pytest.approx(energy_direct)
+ assert_allclose(gradient, gradient_direct)
+ assert obj.last_pair[0][0] == pytest.approx(-1.3)
+ assert obj.last_pair[1][0] == pytest.approx(-1.1)
+ assert obj.last_energy == pytest.approx(energy_direct)
+
+
+def test_mecp_objective_calc_new_honours_geometric_dict_contract():
+ class _FakeScanner:
+ def __call__(self, coords):
+ return ((-1.3, np.ones((3, 3))), (-1.1, np.zeros((3, 3))))
+
+ obj = MECPObjective(_FakeScanner())
+ result = obj.calc_new(np.zeros(9))
+ assert set(result) == {"energy", "gradient"}
+ assert isinstance(result["energy"], float)
+ assert result["gradient"].shape == (9,)
+ assert result["gradient"].ndim == 1
+
+
+def test_mecp_objective_accepts_custom_penalty():
+ calls = []
+
+ def custom(e0, g0, e1, g1, *, sigma, alpha):
+ calls.append((e0, e1, sigma, alpha))
+ return e0 + e1, g0 + g1
+
+ class _FakeScanner:
+ def __call__(self, coords):
+ return ((-1.0, np.zeros(2)), (-2.0, np.ones(2)))
+
+ obj = MECPObjective(_FakeScanner(), sigma=1.0, alpha=0.0, penalty=custom)
+ e, g = obj(None)
+ assert e == pytest.approx(-3.0)
+ assert_allclose(g, np.ones(2))
+ assert calls and calls[0] == (-1.0, -2.0, 1.0, 0.0)
+
+
+# ---------------------------------------------------------------------------
+# Paired-engine mechanics (one SCF + one ADC, two surfaces).
+# ---------------------------------------------------------------------------
+
+def _assert_independent_match(scanner):
+ states = adcc.run_adc(scanner.last_scf,
+ **scanner.paired_target.kwargs())
+ return states
+
+
+def test_paired_excited_scanner_returns_distinct_energies_and_gradients():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="index", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+
+ (e_lo, g_lo), (e_hi, g_hi) = scanner(scanner.initial_coords)
+
+ assert e_lo <= e_hi
+ assert g_lo.shape == (3, 3)
+ assert g_hi.shape == (3, 3)
+ assert np.all(np.isfinite(g_lo)) and np.all(np.isfinite(g_hi))
+
+ # Each channel matches an independent nuclear_gradient eval on the selected
+ # excitation. The paired call selected roots (0, 1) (follow=="index" first
+ # call); energies/gradients must agree elementwise.
+ states = adcc.run_adc(scanner.last_scf, method="adc2", n_singlets=3,
+ conv_tol=1e-9)
+ energies_lo = states.excitations[0].total_energy
+ energies_hi = states.excitations[1].total_energy
+ assert e_lo == pytest.approx(energies_lo, abs=1e-7)
+ assert e_hi == pytest.approx(energies_hi, abs=1e-7)
+ grad_lo = adcc.nuclear_gradient(states.excitations[0],
+ eri_contraction="full_ao").total
+ grad_hi = adcc.nuclear_gradient(states.excitations[1],
+ eri_contraction="full_ao").total
+ assert_allclose(g_lo, grad_lo, atol=1e-7)
+ assert_allclose(g_hi, grad_hi, atol=1e-7)
+
+
+def test_paired_calc_new_returns_both_surfaces():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="index", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+
+ result = scanner.calc_new(scanner.initial_coords.ravel())
+ assert set(result) >= {"energies", "gradients", "energy_lower",
+ "energy_upper", "gradient_lower", "gradient_upper"}
+ assert result["energies"].shape == (2,)
+ assert result["gradients"].shape == (2, 3, 3)
+ assert result["energy_lower"] <= result["energy_upper"]
+ assert result["gradient_lower"].shape == (9,)
+ assert result["gradient_upper"].shape == (9,)
+
+
+def test_paired_excited_overlap_tracking_seeds_and_tracks_both_slots():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="overlap", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+
+ # First call seeds both slots from positional indices; no tracking yet.
+ scanner(scanner.initial_coords)
+ assert scanner.last_trackings == (None, None)
+ assert scanner.previous_indices == (0, 1)
+
+ # Second call at the same geometry: both slots track via overlap.
+ scanner(scanner.initial_coords)
+ slot0, slot1 = scanner.last_trackings
+ assert slot0 is not None and slot1 is not None
+ assert slot0.index != slot1.index # distinctness guard
+ assert slot0.best_score == pytest.approx(1.0, abs=1e-3)
+ assert slot1.best_score == pytest.approx(1.0, abs=1e-3)
+ assert scanner.previous_indices == (slot0.index, slot1.index)
+
+
+def test_paired_mecp_scanner_returns_ground_and_excited_surfaces():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower="mp2", upper=0, n_singlets=3,
+ follow="index", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+
+ (e_lo, g_lo), (e_hi, g_hi) = scanner(scanner.initial_coords)
+
+ assert e_lo <= e_hi
+ assert g_lo.shape == (3, 3) and g_hi.shape == (3, 3)
+
+ # Lower channel is the MP2 ground state: matches LazyMp + nuclear_gradient.
+ mp = adcc.LazyMp(adcc.ReferenceState(scanner.last_scf))
+ assert e_lo == pytest.approx(mp.energy(2), abs=1e-7)
+ exp = adcc.nuclear_gradient(mp, eri_contraction="full_ao").total
+ assert_allclose(g_lo, exp, atol=1e-7)
+
+ # Upper channel is the excited root: only the excited side is tracked.
+ assert scanner.last_excitations[0] is None
+ assert scanner.last_excitations[1] is not None
+ assert scanner.last_trackings == (None, None) # index mode, first call
+
+
+# ---------------------------------------------------------------------------
+# Distinctness guard with synthetic near-degenerate candidates.
+# ---------------------------------------------------------------------------
+
+class _FakeAoOperator:
+ def __init__(self, matrix):
+ self._matrix = matrix
+
+ def to_ndarray(self):
+ return self._matrix
+
+
+class _FakeExcitation:
+ def __init__(self, index, transition_dm, state_diffdm, omega):
+ self.index = index
+ self.transition_dm_ao = _FakeAoOperator(transition_dm)
+ self.state_diffdm_ao = _FakeAoOperator(state_diffdm)
+ self.excitation_energy = omega
+
+
+def test_distinctness_guard_keeps_slots_apart_under_collapse():
+ # Both previous descriptors would, independently, best-match candidate 0.
+ # The joint selection must keep the two roots distinct: one slot takes 0,
+ # the other takes the next-best candidate.
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 0), # same seed is fine here
+ follow="overlap", conv_tol=1e-9,
+ )
+ nao = scanner.base_mol.nao
+
+ def projector(i):
+ m = np.zeros((nao, nao))
+ m[i, i] = 1.0
+ return m
+
+ # Both slots previously followed the *same* character (a collapse scenario).
+ seed_transition = projector(0)
+ seed_diffdm = projector(1)
+ shared_descriptor = _TrackingDescriptor(
+ mol=scanner.base_mol.copy(),
+ transition_dm=seed_transition, state_diffdm=seed_diffdm,
+ )
+
+ # Three candidates: index 0 carries the seeded character, indices 1 and 2
+ # carry orthogonal characters. The guard must not let both slots pick 0.
+ candidates = [
+ _FakeExcitation(0, seed_transition.copy(), seed_diffdm.copy(), 0.20),
+ _FakeExcitation(1, projector(2), projector(3), 0.21),
+ _FakeExcitation(2, projector(4), projector(5), 0.22),
+ ]
+ states = types.SimpleNamespace(excitations=candidates)
+
+ scanner.previous_descriptors = (shared_descriptor, shared_descriptor)
+ scanner.previous_indices = (0, 1)
+
+ chosen, _descriptors, _trackings = scanner._select_excitations(
+ states, scanner.base_mol.copy()
+ )
+
+ assert len(set(chosen)) == 2 # distinct
+ assert 0 in chosen # the best-matching candidate is taken once
+ assert scanner.last_trackings is not None
+
+
+def test_distinctness_guard_raises_when_too_few_states():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), follow="index",
+ )
+ states = types.SimpleNamespace(excitations=[object()]) # only one root
+ with pytest.raises(RuntimeError, match="at least two excited states"):
+ scanner._select_excitations(states, scanner.base_mol.copy())
+
+
+def test_paired_scanner_no_excitations_raises():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), follow="index",
+ )
+ states = types.SimpleNamespace(excitations=[])
+ with pytest.raises(RuntimeError, match="did not return any excited states"):
+ scanner._select_excitations(states, scanner.base_mol.copy())
+
+
+# ---------------------------------------------------------------------------
+# Paired-target normalisation / validation.
+# ---------------------------------------------------------------------------
+
+def test_paired_target_dataclasses_validate_index_arity():
+ PairedExcitedStateTarget(method="adc2", state_indices=(0, 1)) # ok
+ with pytest.raises(ValueError, match="exactly two"):
+ PairedExcitedStateTarget(method="adc2", state_indices=(0, 1, 2))
+
+
+def test_paired_ground_target_enforces_mp2():
+ PairedGroundExcitedStateTarget(method="adc2", state_index=0) # ok, level=2
+ with pytest.raises(NotImplementedError, match="MP2 ground-state"):
+ PairedGroundExcitedStateTarget(method="adc2", state_index=0, level=3)
+
+
+def test_paired_mecp_lower_mp3_raises():
+ with pytest.raises(NotImplementedError, match="MP2 ground-state"):
+ adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower="mp3", upper=0,
+ )
+
+
+def test_paired_mecp_upper_must_be_int():
+ with pytest.raises(TypeError, match="upper must be an int"):
+ adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower="mp2", upper="0",
+ )
+
+
+def test_paired_meci_states_wrong_length_raises():
+ with pytest.raises(ValueError, match="exactly two"):
+ adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1, 2),
+ )
+
+
+def test_paired_meci_missing_method_raises():
+ with pytest.raises(ValueError, match="ADC method is required"):
+ adcc.PairedStateGradientScanner(_h2o_scf(), states=(0, 1))
+
+
+def test_paired_mecp_requires_method():
+ with pytest.raises(ValueError, match="ADC method is required"):
+ adcc.PairedStateGradientScanner(_h2o_scf(), lower="mp2", upper=0)
+
+
+def test_paired_missing_target_raises():
+ with pytest.raises(ValueError, match="paired target"):
+ adcc.PairedStateGradientScanner(_h2o_scf())
+
+
+def test_paired_run_adc_kwargs_forwarded_with_native_names():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=2, conv_tol=1e-7,
+ follow="index",
+ )
+ kwargs = scanner.paired_target.kwargs()
+ assert kwargs["method"] == "adc2"
+ assert kwargs["n_singlets"] == 2
+ assert kwargs["conv_tol"] == pytest.approx(1e-7)
+ assert "output" not in kwargs
+
+
+def test_paired_state_index_out_of_range_raises():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 5), n_singlets=2,
+ follow="index", conv_tol=1e-9,
+ )
+ with pytest.raises(ValueError, match="out of range"):
+ scanner(scanner.initial_coords)
+
+
+def test_paired_invalid_follow_raises_eagerly():
+ with pytest.raises(ValueError, match="overlap.*index"):
+ adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), follow="bogus",
+ )
+
+
+def test_paired_validates_coordinate_shape():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="index", conv_tol=1e-9,
+ )
+ with pytest.raises(ValueError, match="Expected coordinates"):
+ scanner(np.zeros((2, 3)))
+
+
+def test_paired_requires_pyscf_scf_object():
+ with pytest.raises(TypeError, match="PySCF SCF object"):
+ adcc.PairedStateGradientScanner(
+ _h2o_scf().mol, method="adc2", states=(0, 1),
+ )
+
+
+def test_paired_unconverged_scf_raises():
+ from pyscf import scf
+ mf = scf.RHF(_h2o_scf().mol)
+ mf.max_cycle = 1
+ mf.conv_tol = 1e-12
+ mf.conv_tol_grad = 1e-10
+ scanner = adcc.PairedStateGradientScanner(
+ mf, method="adc2", states=(0, 1), n_singlets=3, follow="index",
+ )
+ with pytest.raises(RuntimeError, match="did not converge"):
+ scanner(scanner.initial_coords)
+
+
+# ---------------------------------------------------------------------------
+# geomeTRIC oracle cross-check (only when geometric is importable).
+#
+# This calls geomeTRIC's *real* :class:`geometric.engine.ConicalIntersection`
+# engine with fake single-point sub-engines returning controlled energies and
+# gradients, then asserts :func:`mecp_penalty` reproduces geomeTRIC's own
+# penalty objective and gradient to machine precision. This is a direct
+# cross-check against the implementation that already exists upstream, not a
+# re-port of its formula.
+# ---------------------------------------------------------------------------
+
+def _geom_conical_intersection(energies, grads, sigma, alpha):
+ """Build a live geomeTRIC ConicalIntersection over fake single-point engines.
+
+ Returns the ``(energy, gradient)`` geomeTRIC computes for the given
+ controlled per-state energies and gradients. The fake sub-engines subclass
+ :class:`geometric.engine.Engine` and return canned single-point dicts from
+ ``calc_new`` so the real penalty machinery in ``ConicalIntersection.calc_new``
+ runs end-to-end.
+ """
+ import tempfile
+ from geometric.engine import ConicalIntersection, Engine
+
+ class _FakeMol:
+ def __len__(self):
+ return 1
+
+ class _FakeEngine(Engine):
+ def __init__(self, molecule, energy, gradient):
+ self._energy = float(energy)
+ self._gradient = np.asarray(gradient, dtype=float)
+ super().__init__(molecule)
+
+ def calc_new(self, coords, dirname):
+ return {"energy": self._energy,
+ "gradient": self._gradient.copy()}
+
+ mol = _FakeMol()
+ engines = [_FakeEngine(mol, e, g) for e, g in zip(energies, grads)]
+ ci = ConicalIntersection(mol, engines, sigma, alpha)
+ coords = np.zeros(np.asarray(grads[0]).size)
+ with tempfile.TemporaryDirectory() as dnm:
+ out = ci.calc(coords, dnm)
+ return out["energy"], np.asarray(out["gradient"])
+
+
+def test_penalty_matches_geometric_oracle():
+ import importlib.util
+ if importlib.util.find_spec("geometric") is None:
+ pytest.skip("geometric not installed")
+ rng = np.random.default_rng(123)
+ energies = [-1.3, -1.08]
+ grads = [rng.normal(size=6), rng.normal(size=6)]
+ for sigma, alpha in [(3.5, 0.025), (5.0, 0.1), (1.0, 0.0)]:
+ e_oracle, g_oracle = _geom_conical_intersection(
+ energies, grads, sigma, alpha,
+ )
+ e_adcc, g_adcc = mecp_penalty(
+ energies[0], grads[0], energies[1], grads[1],
+ sigma=sigma, alpha=alpha,
+ )
+ assert e_adcc == pytest.approx(e_oracle, abs=1e-13)
+ assert_allclose(g_adcc, g_oracle, atol=1e-13)
+
+
+def test_penalty_matches_geometric_oracle_at_degeneracy():
+ # At exact degeneracy the penalty and its gradient vanish; the oracle must
+ # reduce to the average surface, matching mecp_penalty exactly.
+ import importlib.util
+ if importlib.util.find_spec("geometric") is None:
+ pytest.skip("geometric not installed")
+ rng = np.random.default_rng(7)
+ energies = [-1.2, -1.2]
+ grads = [rng.normal(size=6), rng.normal(size=6)]
+ e_oracle, g_oracle = _geom_conical_intersection(energies, grads, 3.5, 0.025)
+ e_adcc, g_adcc = mecp_penalty(
+ energies[0], grads[0], energies[1], grads[1], sigma=3.5, alpha=0.025,
+ )
+ assert e_adcc == pytest.approx(e_oracle, abs=1e-13)
+ assert_allclose(g_adcc, g_oracle, atol=1e-13)
+ # Sanity: the oracle itself collapsed to the average at the seam.
+ assert e_oracle == pytest.approx(-1.2, abs=1e-13)
+ assert_allclose(g_oracle, 0.5 * (np.asarray(grads[0]) + np.asarray(grads[1])))
diff --git a/docs/gradients.rst b/docs/gradients.rst
index 6acb5924..5dac98cb 100644
--- a/docs/gradients.rst
+++ b/docs/gradients.rst
@@ -90,9 +90,64 @@ The same object also implements geomeTRIC's custom-engine ``calc_new`` protocol:
# result == {"energy": energy, "gradient": gradient.ravel()}
geomeTRIC remains an optional dependency; the plain scanner only requires the
-PySCF backend. Minimum-energy crossing point workflows can be built from two
-scanner targets because geomeTRIC's penalty-constrained formulation only needs
-the two state energies and gradients, not derivative couplings.
+PySCF backend.
+
+.. _gradients-mecp:
+
+Minimum-energy crossing points (MECP/MECI)
+------------------------------------------
+
+The :class:`adcc.PairedStateGradientScanner` evaluates *two* electronic
+surfaces at one geometry from a single SCF and a single ADC (or, for a
+ground/excited MECP, one ``LazyMp`` plus one ADC) and returns both
+``(energy, gradient)`` pairs. It reuses the single-surface scanner's SCF
+lifecycle and AO density-overlap root tracking, and adds a **joint distinct root
+selection** with a distinctness guard so the two followed roots never collapse
+onto the same candidate state -- the key complication when two surfaces become
+degenerate at a conical-intersection seam.
+
+No derivative (non-adiabatic) couplings are required: the penalty-function
+family of MECP/MECI optimisers only needs the two state energies and their
+gradients. Two normalised target forms are supported. For an excited/excited
+MECI pair the two roots come from one ADC solve via the ``states=(i, j)``
+convenience::
+
+ paired = adcc.PairedStateGradientScanner(
+ scfres,
+ method="adc2",
+ states=(0, 1), # two tracked excited states (MECI)
+ n_singlets=5,
+ follow="overlap",
+ )
+ (e_lower, g_lower), (e_upper, g_upper) = paired(mol.atom_coords())
+
+For a ground/excited MECP pair the ground surface is always MP2 (the only level
+with an analytic gradient in adcc and is evaluated via :class:`adcc.LazyMp`)
+paired with a single tracked excited root::
+
+ paired = adcc.PairedStateGradientScanner(
+ scfres,
+ method="adc2",
+ lower="mp2", upper=0, # MP2 ground + excited root 0 (MECP)
+ n_singlets=5,
+ follow="overlap",
+ )
+
+To drive a geomeTRIC optimisation, wrap the paired scanner in an
+:class:`adcc.MECPObjective`. This combines the two surfaces into a single
+penalty ``(energy, gradient)`` using the smoothed Levine--Coe--Martinez penalty
+(the same form as geomeTRIC's built-in conical-intersection engine), with the
+``sigma`` and ``alpha`` parameters tuning the penalty strength and seam
+smoothing::
+
+ from pyscf.geomopt import as_pyscf_method, geometric_solver
+
+ objective = adcc.MECPObjective(paired)
+ method = as_pyscf_method(mol, lambda m: objective(m.atom_coords(unit="Bohr")))
+ mol_ci = geometric_solver.optimize(method, maxsteps=20)
+
+A complete worked example is provided in
+``examples/optimization/pyscf_adcc_mecp.py``.
The two-electron term can be evaluated with two different strategies, selected
through the ``eri_contraction`` keyword (see
diff --git a/examples/optimization/pyscf_adcc_mecp.py b/examples/optimization/pyscf_adcc_mecp.py
new file mode 100644
index 00000000..1705d2a1
--- /dev/null
+++ b/examples/optimization/pyscf_adcc_mecp.py
@@ -0,0 +1,101 @@
+#!/usr/bin/env python3
+"""Minimum-energy conical intersection (MECI) optimisation with adcc gradients.
+
+This example locates a conical-intersection seam point of twisted ethylene
+(C2H4), the textbook minimal MECI system, using the paired-state scanner and the
+Levine--Coe--Martinez penalty objective. The
+:class:`adcc.PairedStateGradientScanner` evaluates two excited-state surfaces at
+one geometry from a single SCF + single
+ADC, and :class:`adcc.MECPObjective` combines them into one
+``(energy, gradient)`` for geomeTRIC. **No derivative couplings** are required.
+
+geomeTRIC is optional; install it separately, e.g. ``pip install geometric`` or
+``pip install adcc[geomopt]`` once the optional extra is available.
+"""
+
+import numpy as np
+
+import adcc
+from pyscf import gto, scf
+from pyscf.geomopt import as_pyscf_method, geometric_solver
+
+
+def twisted_ethylene(twist_deg=80.0, basis="sto-3g"):
+ """Build a twisted ethylene PySCF SCF near the S1/S0 crossing seam."""
+ tw = np.deg2rad(twist_deg)
+ c1 = np.array([0.0, 0.0, 0.0])
+ c2 = np.array([1.34, 0.0, 0.0])
+ h1l = c1 + np.array([0.0, 0.63, 0.0])
+ h2l = c1 + np.array([0.0, -0.63, 0.0])
+
+ def twist(p, angle, pivot):
+ r = p - pivot
+ rot = np.array([[1.0, 0.0, 0.0],
+ [0.0, np.cos(angle), -np.sin(angle)],
+ [0.0, np.sin(angle), np.cos(angle)]])
+ return pivot + rot @ r
+
+ h1r = twist(c2 + np.array([0.0, 0.63, 0.0]), tw, c2)
+ h2r = twist(c2 + np.array([0.0, -0.63, 0.0]), tw, c2)
+ atoms = ["C", "C", "H", "H", "H", "H"]
+ coords = np.stack([c1, c2, h1l, h2l, h1r, h2r])
+ atom_str = "\n".join(
+ f"{sym} {xyz[0]:.8f} {xyz[1]:.8f} {xyz[2]:.8f}"
+ for sym, xyz in zip(atoms, coords)
+ )
+ mol = gto.M(
+ atom=atom_str, basis=basis, unit="Angstrom",
+ symmetry=False, verbose=0, parse_arg=False,
+ )
+ mf = scf.RHF(mol)
+ mf.conv_tol = 1e-10
+ mf.conv_tol_grad = 1e-7
+ return mf
+
+
+def print_geometry(title, mol):
+ print(f"\n{title} (Angstrom):")
+ for i, xyz in enumerate(mol.atom_coords(unit="Angstrom")):
+ print(f"{mol.atom_symbol(i):2s} "
+ f"{xyz[0]:16.10f} {xyz[1]:16.10f} {xyz[2]:16.10f}")
+
+
+if __name__ == "__main__":
+ scfres = twisted_ethylene(twist_deg=80.0)
+ mol = scfres.mol
+ print_geometry("Initial twisted geometry", mol)
+
+ paired = adcc.PairedStateGradientScanner(
+ scfres,
+ method="adc2",
+ states=(0, 1), # two tracked excited states (MECI)
+ n_singlets=4,
+ follow="overlap",
+ conv_tol=1e-8,
+ gradient_kwargs={"eri_contraction": "full_ao", "conv_tol": 1e-7},
+ )
+ # The default penalty uses the smoothed Levine--Coe--Martinez form, the same
+ # formulation as geomeTRIC's built-in conical-intersection engine.
+ objective = adcc.MECPObjective(paired)
+
+ def energy_and_gradient(mol_at_step):
+ energy, gradient = objective(mol_at_step.atom_coords(unit="Bohr"))
+ e_lo, e_hi = objective.last_pair[0][0], objective.last_pair[1][0]
+ gnorm = float((gradient ** 2).sum() ** 0.5)
+ print(f"E_pen = {energy:.12f} Eh, |g| = {gnorm:.6e} Eh/Bohr, "
+ f"gap = {abs(e_hi - e_lo):.6e} Eh")
+ return energy, gradient
+
+ method = as_pyscf_method(mol, energy_and_gradient)
+ mol_ci = geometric_solver.optimize(
+ method,
+ maxsteps=40,
+ convergence_grms=1e-4,
+ convergence_gmax=2e-4,
+ )
+ print_geometry("MECI geometry", mol_ci)
+
+ # Report the final two surfaces at the located geometry.
+ (e_lo, _), (e_hi, _) = paired(mol_ci.atom_coords(unit="Bohr"))
+ print(f"\nFinal surface energies: lower = {e_lo:.10f} Eh, "
+ f"upper = {e_hi:.10f} Eh, gap = {abs(e_hi - e_lo):.2e} Eh")
From b417787e69c4268c394e4b0ff9039f301a85ecbb Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Mon, 22 Jun 2026 14:48:38 +0200
Subject: [PATCH 3/9] Harden MECP/MECI penalty scanner and convert example to
S0/S1 MECP
- mecp_penalty: short-circuit the alpha=0 / exact-degeneracy singular point to
the average surface (was ZeroDivisionError); guard against non-finite
objective output; scope the docstring to alpha > 0.
- paired_scanner: validate np.isfinite energies/gradients before the energy
argsort (NaN was silently reordering the lower/upper slots); make
tracking_min_gap raise (not silently inert) when only two excited states are
computed on both the MECI joint and MECP single-root paths; correct the
wrong inline comment about an always-present alternative.
- geomopt_mecp_test: add 22 fast unit tests for the singular degeneracy guard,
non-finite rejection, paired/MECP threshold enforcement, two-state min_gap
behaviour, MECP excited-root overlap tracking, per-channel diagnostics,
energy-sort reorder, MECP lower normalisation branches and out-of-range
upper, index-mode collapse, and the joint seam-bias tie-break.
- pyscf_adcc_mecp example: locate the S0/S1 MECP (MP2 ground + first tracked
excited singlet) instead of the S1/S2 MECI.
---
adcc/gradients/mecp.py | 36 ++-
adcc/gradients/paired_scanner.py | 34 +-
adcc/tests/geomopt_mecp_test.py | 379 +++++++++++++++++++++++
examples/optimization/pyscf_adcc_mecp.py | 39 +--
4 files changed, 458 insertions(+), 30 deletions(-)
diff --git a/adcc/gradients/mecp.py b/adcc/gradients/mecp.py
index bfdb7659..9956e281 100644
--- a/adcc/gradients/mecp.py
+++ b/adcc/gradients/mecp.py
@@ -88,7 +88,11 @@ def mecp_penalty(e_lower, g_lower, e_upper, g_upper, *, sigma=DEFAULT_SIGMA,
At exact degeneracy (``e_upper == e_lower``) the penalty and its gradient
vanish, leaving ``E_obj = E_avg`` and ``G_obj = G_avg`` -- the objective
- reduces to an average-surface optimisation pinned to the crossing.
+ reduces to an average-surface optimisation pinned to the crossing. This
+ includes the ``alpha == 0`` raw energy-difference mode: the limiting
+ penalty ``sigma * E_dif`` tends to zero as ``E_dif -> 0`` even though the
+ un-guarded expression ``sigma * E_dif**2 / (E_dif + 0)`` is ``0/0`` there,
+ so the exact-degeneracy point is short-circuited explicitly below.
Parameters
----------
@@ -123,17 +127,33 @@ def mecp_penalty(e_lower, g_lower, e_upper, g_upper, *, sigma=DEFAULT_SIGMA,
g_avg = 0.5 * (g_lower + g_upper)
# Smoothed penalty (Levine--Coe--Martinez / geomeTRIC ConicalIntersection).
- # At e_dif == 0 the penalty term and its gradient vanish by construction.
- denom = e_dif + alpha
- e_pen = sigma * (e_dif * e_dif) / (denom * _N_STATES2)
- g_pen_scale = (
- sigma * (e_dif * e_dif + 2.0 * alpha * e_dif)
- / (denom * denom * _N_STATES2)
- )
+ # At exact degeneracy the penalty term and its gradient vanish by
+ # construction (the LCM objective is continuously extended to the seam).
+ # This is the only singular point: with alpha == 0 the raw
+ # energy-difference mode would otherwise evaluate 0.0 / .0 here, so we
+ # short-circuit the seam explicitly. The guard also covers the (only
+ # theoretical) sub-DBL_MIN underflow where e_dif**2 flushes to zero while
+ # alpha == 0.
+ if e_dif == 0.0 or (alpha == 0.0 and abs(e_dif) < 1e-300):
+ e_pen = 0.0
+ g_pen_scale = 0.0
+ else:
+ denom = e_dif + alpha
+ e_pen = sigma * (e_dif * e_dif) / (denom * _N_STATES2)
+ g_pen_scale = (
+ sigma * (e_dif * e_dif + 2.0 * alpha * e_dif)
+ / (denom * denom * _N_STATES2)
+ )
g_pen = g_pen_scale * g_dif
e_obj = e_avg + e_pen
g_obj = g_avg + g_pen
+ if not np.isfinite(e_obj) or not np.all(np.isfinite(g_obj)):
+ raise RuntimeError(
+ "mecp_penalty produced a non-finite objective (E_obj="
+ f"{e_obj!r}, e_lower={e_lower!r}, e_upper={e_upper!r}); a "
+ "non-finite state energy or gradient reached the penalty."
+ )
return float(e_obj), np.asarray(g_obj, dtype=float)
diff --git a/adcc/gradients/paired_scanner.py b/adcc/gradients/paired_scanner.py
index 611ee42a..8c81413c 100644
--- a/adcc/gradients/paired_scanner.py
+++ b/adcc/gradients/paired_scanner.py
@@ -243,6 +243,17 @@ def __call__(self, coords):
self.last_gradients = tuple(grads)
self.last_energies = tuple(energies)
+ # Guard against non-finite surface results *before* the energy sort:
+ # np.argsort silently reorders slots on NaN, which would otherwise let a
+ # borderline ADC solve hand misleading (lower, upper) pairs -- and a
+ # NaN objective -- to geomeTRIC.
+ for slot, e in enumerate(energies):
+ if not np.isfinite(e) or not np.all(np.isfinite(grads[slot])):
+ raise RuntimeError(
+ f"Non-finite surface result on paired slot {slot} "
+ f"(energy={e!r}); cannot drive a penalty optimisation."
+ )
+
order = np.argsort(energies)
lo, hi = int(order[0]), int(order[1])
self.last_pair_order = (lo, hi)
@@ -504,12 +515,20 @@ def _select_excitations(self, states, mol):
continue
other = chosen[1 - slot]
# Gap to the best *alternative* root available to this slot, i.e.
- # excluding the partner's selected root (the distinctness guard
- # means there is always at least one other candidate here).
+ # excluding the partner's selected root. When exactly two excited
+ # states were computed (n == 2) each slot is forced onto the
+ # partner's leftover root, so no per-channel "best vs second-best"
+ # gap is measurable here: ``remaining`` is empty.
remaining = [k for k in range(n) if k != idx and k != other]
if remaining:
second = max(remaining, key=lambda k: sv[k])
gap = float(sv[idx] - sv[second])
+ elif self.tracking_min_gap > 0:
+ raise RuntimeError(
+ f"tracking_min_gap > 0 needs at least three computed "
+ f"excited states to measure a per-channel gap for paired "
+ f"slot {slot} (only {n} available)."
+ )
else:
gap = float("inf")
prev_idx = self.previous_indices[slot]
@@ -568,8 +587,15 @@ def _select_excited_root(self, states, mol, prev_descriptor, prev_index,
])
order = np.argsort(scores)[::-1]
best = int(order[0])
- gap = (float(scores[best] - scores[int(order[1])])
- if len(order) > 1 else float("inf"))
+ if len(order) > 1:
+ gap = float(scores[best] - scores[int(order[1])])
+ elif self.tracking_min_gap > 0:
+ raise RuntimeError(
+ "tracking_min_gap > 0 needs at least two computed excited "
+ f"states to measure a gap (only {n} available)."
+ )
+ else:
+ gap = float("inf")
switched = prev_index is not None and best != prev_index
tracking = TrackingResult(
index=best, scores=scores, best_score=float(scores[best]), gap=gap,
diff --git a/adcc/tests/geomopt_mecp_test.py b/adcc/tests/geomopt_mecp_test.py
index 0ee707de..ea9f1e72 100644
--- a/adcc/tests/geomopt_mecp_test.py
+++ b/adcc/tests/geomopt_mecp_test.py
@@ -612,3 +612,382 @@ def test_penalty_matches_geometric_oracle_at_degeneracy():
# Sanity: the oracle itself collapsed to the average at the seam.
assert e_oracle == pytest.approx(-1.2, abs=1e-13)
assert_allclose(g_oracle, 0.5 * (np.asarray(grads[0]) + np.asarray(grads[1])))
+
+
+# ---------------------------------------------------------------------------
+# Exact-degeneracy guard at alpha == 0 (findings 1 & 5) and finiteness guards.
+# ---------------------------------------------------------------------------
+
+def test_penalty_alpha_zero_exact_degeneracy_collapses_to_average():
+ # At (alpha == 0, e_dif == 0) the raw energy-difference mode would
+ # otherwise evaluate 0.0 / 0.0; the guard short-circuits to the average
+ # surface so the penalty vanishes exactly, matching the docstring.
+ g_lo = np.array([[1.0, 0.0, 0.0]])
+ g_hi = np.array([[0.0, 2.0, 0.0]])
+ e = -1.2
+ energy, gradient = mecp_penalty(e, g_lo, e, g_hi, sigma=2.0, alpha=0.0)
+ assert energy == pytest.approx(e)
+ assert_allclose(gradient, 0.5 * (g_lo + g_hi))
+
+
+def test_penalty_alpha_zero_near_degeneracy_is_finite():
+ # The sub-DBL_MIN underflow guard (-- abs(e_dif) < 1e-300 with alpha == 0)
+ # must also short-circuit rather than produce 0.0 / 0.0.
+ tiny = 1e-310 # e_dif ** 2 flushes to zero in float64
+ energy, gradient = mecp_penalty(
+ -1.2, np.zeros(3), -1.2 + tiny, np.zeros(3), sigma=2.0, alpha=0.0,
+ )
+ assert np.isfinite(energy)
+
+ def efun(off):
+ return mecp_penalty(-1.2, np.zeros(3), -1.2 + off, np.zeros(3),
+ sigma=2.0, alpha=0.0)[0]
+
+ # The limit from outside the guard (alpha == 0) is the average as off -> 0.
+ assert efun(tiny) == pytest.approx(-1.2, abs=1e-6)
+
+
+def test_penalty_rejects_non_finite_energies():
+ g = np.zeros(3)
+ with pytest.raises(RuntimeError, match="non-finite"):
+ mecp_penalty(float("nan"), g, -1.0, g)
+ with pytest.raises(RuntimeError, match="non-finite"):
+ mecp_penalty(-1.0, g, float("inf"), g)
+
+
+def test_penalty_rejects_non_finite_gradients():
+ with pytest.raises(RuntimeError, match="non-finite"):
+ mecp_penalty(-1.3, np.array([float("nan"), 0.0, 0.0]), -1.1,
+ np.zeros(3))
+
+
+# ---------------------------------------------------------------------------
+# Paired threshold enforcement (finding 3) for MECI and MECP overlap tracking.
+# ---------------------------------------------------------------------------
+
+def test_paired_overlap_tracking_below_min_score_raises():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="overlap", tracking_min_score=2.0, conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ scanner(scanner.initial_coords) # seed
+ with pytest.raises(RuntimeError, match="below threshold"):
+ scanner(scanner.initial_coords)
+
+
+def test_paired_overlap_tracking_ambiguous_gap_raises():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="overlap", tracking_min_gap=10.0, conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ scanner(scanner.initial_coords) # seed
+ with pytest.raises(RuntimeError, match="ambiguous"):
+ scanner(scanner.initial_coords)
+
+
+def test_mecp_overlap_tracking_below_min_score_raises():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower="mp2", upper=0, n_singlets=3,
+ follow="overlap", tracking_min_score=2.0, conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ scanner(scanner.initial_coords) # seed (ground never tracked)
+ with pytest.raises(RuntimeError, match="below threshold"):
+ scanner(scanner.initial_coords)
+
+
+def test_mecp_overlap_tracking_ambiguous_gap_raises():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower="mp2", upper=0, n_singlets=3,
+ follow="overlap", tracking_min_gap=10.0, conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ scanner(scanner.initial_coords) # seed
+ with pytest.raises(RuntimeError, match="ambiguous"):
+ scanner(scanner.initial_coords)
+
+
+# ---------------------------------------------------------------------------
+# tracking_min_gap with only two states (finding 4).
+# ---------------------------------------------------------------------------
+
+def test_tracking_min_gap_raises_when_only_two_excited_states():
+ # n_singlets == 2 leaves no per-channel alternative to compare a gap
+ # against; requesting ambiguity detection must surface that clearly.
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=2,
+ follow="overlap", tracking_min_gap=0.5, conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ scanner(scanner.initial_coords) # seed (no tracking yet)
+ with pytest.raises(RuntimeError, match="at least three computed excited"):
+ scanner(scanner.initial_coords)
+
+
+def test_tracking_min_gap_inert_at_two_states_with_default():
+ # With the default tracking_min_gap == 0 the two-state case does not raise;
+ # the per-channel gap is reported as inf (no measurable alternative).
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=2,
+ follow="overlap", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ scanner(scanner.initial_coords)
+ scanner(scanner.initial_coords)
+ slot0, slot1 = scanner.last_trackings
+ assert slot0 is not None and slot1 is not None
+ assert slot0.gap == float("inf")
+ assert slot1.gap == float("inf")
+ assert slot0.index != slot1.index
+
+
+# ---------------------------------------------------------------------------
+# MECP excited-root overlap tracking (finding 2) and per-channel diagnostics.
+# ---------------------------------------------------------------------------
+
+def test_mecp_overlap_tracking_seeds_and_tracks_excited_root():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower="mp2", upper=0, n_singlets=3,
+ follow="overlap", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ # Seed call: ground is never tracked; the excited side falls back to its
+ # positional seed index because there is no previous descriptor yet.
+ scanner(scanner.initial_coords)
+ assert scanner.last_trackings == (None, None)
+ assert scanner.last_excitations[0] is None
+ assert scanner.last_excitations[1] is not None
+
+ # Second call at the same geometry: the excited side tracks via overlap.
+ scanner(scanner.initial_coords)
+ ground_tr, excited_tr = scanner.last_trackings
+ assert ground_tr is None # ground remains untracked
+ assert excited_tr is not None
+ assert excited_tr.index == 0 # the seeded S1 character
+ assert excited_tr.best_score == pytest.approx(1.0, abs=1e-3)
+ assert excited_tr.previous_index == 0
+ assert excited_tr.switched is False
+
+
+def test_paired_overlap_tracking_records_per_channel_diagnostics():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="overlap", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ scanner(scanner.initial_coords) # seed
+ scanner(scanner.initial_coords) # track both slots via overlap
+ slot0, slot1 = scanner.last_trackings
+ assert slot0 is not None and slot1 is not None
+ assert slot0.index != slot1.index
+ # Per-channel diagnostics that catch a silent root flip.
+ for slot in (slot0, slot1):
+ assert np.isfinite(slot.gap)
+ assert slot.gap >= 0.0
+ assert slot0.previous_index == 0
+ assert slot1.previous_index == 1
+ assert slot0.switched is False
+ assert slot1.switched is False
+
+
+# ---------------------------------------------------------------------------
+# Energy-sorting reorder branch (finding 8) via synthetic surface injection.
+# ---------------------------------------------------------------------------
+
+class _FakeSurface:
+ """A non-LazyMp target with a scalar total_energy and (ignored) extras."""
+
+ def __init__(self, total_energy):
+ self.total_energy = total_energy
+
+
+class _FakeGradResult:
+ def __init__(self, total):
+ self.total = np.asarray(total, dtype=float)
+
+
+def _patched_paired_scanner(monkeypatch, energies):
+ """Build a paired scanner whose __call__ uses fake, injected surfaces.
+
+ ``_run_scf`` and ``_build_target`` are stubbed so the energy path runs on
+ controlled ``total_energy`` values, and the gradient eval is stubbed to a
+ zero gradient. This isolates the energy-sorting / finiteness logic of
+ ``PairedStateGradientScanner.__call__`` from any SCF/ADC machinery.
+ """
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="index", conv_tol=1e-9,
+ )
+ monkeypatch.setattr(scanner, "_run_scf", lambda coords: object())
+ targets = tuple(_FakeSurface(e) for e in energies)
+ monkeypatch.setattr(scanner, "_build_target", lambda scfres: targets)
+ import adcc.gradients as gradients_mod
+ monkeypatch.setattr(
+ gradients_mod, "nuclear_gradient",
+ lambda t, **kwargs: _FakeGradResult(np.zeros((3, 3))),
+ )
+ return scanner
+
+
+def test_paired_call_energy_sorts_and_records_last_pair_order(monkeypatch):
+ # Slot 1 carries the lower energy; __call__ must return it as the lower
+ # surface while keeping last_energies in slot order.
+ e_slot0, e_slot1 = -1.1, -1.3 # slot 1 is lower
+ scanner = _patched_paired_scanner(monkeypatch, (e_slot0, e_slot1))
+ (e_lo, g_lo), (e_hi, g_hi) = scanner(scanner.initial_coords)
+ assert e_lo == pytest.approx(e_slot1)
+ assert e_hi == pytest.approx(e_slot0)
+ assert scanner.last_pair_order == (1, 0) # energy-sorted insertion
+ assert scanner.last_energies == (e_slot0, e_slot1) # slot order preserved
+
+
+def test_paired_call_default_slot_order_when_already_sorted(monkeypatch):
+ scanner = _patched_paired_scanner(monkeypatch, (-1.3, -1.1)) # slot0 lower
+ (e_lo, _), (e_hi, _) = scanner(scanner.initial_coords)
+ assert e_lo == pytest.approx(-1.3)
+ assert e_hi == pytest.approx(-1.1)
+ assert scanner.last_pair_order == (0, 1)
+
+
+# ---------------------------------------------------------------------------
+# Non-finite surface guard in the scanner (finding 7, paired-side).
+# ---------------------------------------------------------------------------
+
+def test_paired_scanner_rejects_non_finite_energy(monkeypatch):
+ import adcc.gradients as gradients_mod
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3, follow="index",
+ )
+ monkeypatch.setattr(scanner, "_run_scf", lambda coords: object())
+ monkeypatch.setattr(
+ scanner, "_build_target",
+ lambda scfres: (_FakeSurface(float("inf")), _FakeSurface(-1.0)),
+ )
+ monkeypatch.setattr(
+ gradients_mod, "nuclear_gradient",
+ lambda t, **kwargs: _FakeGradResult(np.zeros((3, 3))),
+ )
+ with pytest.raises(RuntimeError, match="Non-finite surface result"):
+ scanner(scanner.initial_coords)
+
+
+def test_paired_scanner_rejects_non_finite_gradient(monkeypatch):
+ import adcc.gradients as gradients_mod
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3, follow="index",
+ )
+ monkeypatch.setattr(scanner, "_run_scf", lambda coords: object())
+ monkeypatch.setattr(
+ scanner, "_build_target",
+ lambda scfres: (_FakeSurface(-1.0), _FakeSurface(-2.0)),
+ )
+ monkeypatch.setattr(
+ gradients_mod, "nuclear_gradient",
+ lambda t, **kwargs: _FakeGradResult(np.full((3, 3), float("nan"))),
+ )
+ with pytest.raises(RuntimeError, match="Non-finite surface result"):
+ scanner(scanner.initial_coords)
+
+
+# ---------------------------------------------------------------------------
+# MECP lower normalisation branches and out-of-range upper (finding 10).
+# ---------------------------------------------------------------------------
+
+def test_paired_mecp_lower_as_ground_state_target():
+ from adcc.gradients import GroundStateTarget
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower=GroundStateTarget(level=2), upper=0,
+ )
+ assert isinstance(scanner.paired_target, PairedGroundExcitedStateTarget)
+ assert scanner.paired_target.level == 2
+
+
+def test_paired_mecp_lower_as_int_level():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower=2, upper=0,
+ )
+ assert isinstance(scanner.paired_target, PairedGroundExcitedStateTarget)
+ assert scanner.paired_target.level == 2
+
+
+def test_paired_mecp_lower_bad_type_raises():
+ with pytest.raises(TypeError, match="lower must be"):
+ adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower=1.5, upper=0,
+ )
+
+
+def test_paired_mecp_upper_out_of_range_raises():
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower="mp2", upper=5, n_singlets=2,
+ follow="index", conv_tol=1e-9,
+ )
+ with pytest.raises(ValueError, match="out of range"):
+ scanner(scanner.initial_coords)
+
+
+# ---------------------------------------------------------------------------
+# Distinctness guard on duplicate positional seeds (finding 11).
+# ---------------------------------------------------------------------------
+
+def test_paired_meci_duplicate_index_mode_raises_distinctness():
+ # In pure index mode with states=(0, 0) both slots want the same root;
+ # the distinctness guard must raise rather than silently collapse.
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 0), n_singlets=3,
+ follow="index", conv_tol=1e-9,
+ )
+ with pytest.raises(RuntimeError, match="collapsed onto state"):
+ scanner(scanner.initial_coords)
+
+
+# ---------------------------------------------------------------------------
+# Joint tie-break toward the smaller excitation-energy gap (finding 13).
+# ---------------------------------------------------------------------------
+
+def test_joint_selection_tie_break_prefers_smaller_excitation_gap(monkeypatch):
+ # Both slots track via overlap with two distinct ordered pairs tied on the
+ # summed overlap score but differing in excitation-energy gap; the joint
+ # selection must prefer the pair closer to the seam (smaller gap) even
+ # though it is neither the first-iterated nor the lowest-index pair.
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=4,
+ follow="overlap", conv_tol=1e-9,
+ )
+ nao = scanner.base_mol.nao
+
+ # Score table keyed by the previous descriptor identity + candidate index:
+ # slot0 likes {0, 2}, slot1 likes {1, 3} -> four ordered pairs tie at 1.0.
+ sentinel0, sentinel1 = object(), object()
+ score_map = {
+ id(sentinel0): [1.0, 0.0, 1.0, 0.0],
+ id(sentinel1): [0.0, 1.0, 0.0, 1.0],
+ }
+ omegas = [0.10, 0.20, 0.11, 0.19]
+ # Paired gap analysis: (0,1)=0.10 (0,3)=0.09 (2,1)=0.09 (2,3)=0.08 smallest.
+
+ # _descriptor returns the candidate index so _tracking_score can map it
+ # back to a controlled score without fabricating density matrices.
+ monkeypatch.setattr(
+ scanner, "_descriptor", lambda excitation, mol: excitation.index,
+ )
+ monkeypatch.setattr(
+ scanner, "_tracking_score",
+ lambda prev, current, mol, unavailable:
+ float(score_map[id(prev)][current]),
+ )
+ # Ensure score_vectors are built (follow == "overlap" AND prev is not None).
+ scanner.previous_descriptors = (sentinel0, sentinel1)
+ scanner.previous_indices = (0, 1)
+
+ candidates = [
+ _FakeExcitation(i, np.zeros((nao, nao)), np.zeros((nao, nao)), omegas[i])
+ for i in range(4)
+ ]
+ states = types.SimpleNamespace(excitations=candidates)
+ chosen, _descriptors, _trackings = scanner._select_excitations(
+ states, scanner.base_mol.copy()
+ )
+ assert tuple(chosen) == (2, 3) # smallest excitation-energy-gap pair
diff --git a/examples/optimization/pyscf_adcc_mecp.py b/examples/optimization/pyscf_adcc_mecp.py
index 1705d2a1..25a732d9 100644
--- a/examples/optimization/pyscf_adcc_mecp.py
+++ b/examples/optimization/pyscf_adcc_mecp.py
@@ -1,12 +1,13 @@
#!/usr/bin/env python3
-"""Minimum-energy conical intersection (MECI) optimisation with adcc gradients.
-
-This example locates a conical-intersection seam point of twisted ethylene
-(C2H4), the textbook minimal MECI system, using the paired-state scanner and the
-Levine--Coe--Martinez penalty objective. The
-:class:`adcc.PairedStateGradientScanner` evaluates two excited-state surfaces at
-one geometry from a single SCF + single
-ADC, and :class:`adcc.MECPObjective` combines them into one
+"""Minimum-energy crossing point (MECP) of the S0/S1 seam with adcc gradients.
+
+This example locates the minimum on the crossing seam between the ground state
+S0 and the first excited singlet S1 of twisted ethylene (C2H4) -- the textbook
+S0/S1 MECP / conical-intersection system -- using the paired-state scanner and
+the Levine--Coe--Martinez penalty objective. The
+:class:`adcc.PairedStateGradientScanner` evaluates the ground-state MP2 surface
+and one tracked ADC excited surface at a single geometry from a single SCF +
+single ADC, and :class:`adcc.MECPObjective` combines them into one
``(energy, gradient)`` for geomeTRIC. **No derivative couplings** are required.
geomeTRIC is optional; install it separately, e.g. ``pip install geometric`` or
@@ -50,6 +51,7 @@ def twist(p, angle, pivot):
mf = scf.RHF(mol)
mf.conv_tol = 1e-10
mf.conv_tol_grad = 1e-7
+ mf.max_cycles = 250
return mf
@@ -61,18 +63,18 @@ def print_geometry(title, mol):
if __name__ == "__main__":
- scfres = twisted_ethylene(twist_deg=80.0)
+ scfres = twisted_ethylene(twist_deg=80.0, basis="6-31g")
mol = scfres.mol
print_geometry("Initial twisted geometry", mol)
paired = adcc.PairedStateGradientScanner(
scfres,
method="adc2",
- states=(0, 1), # two tracked excited states (MECI)
- n_singlets=4,
+ lower="mp2", upper=0, # S0 (MP2 ground) + S1 (1st tracked excited)
+ n_singlets=3,
follow="overlap",
- conv_tol=1e-8,
- gradient_kwargs={"eri_contraction": "full_ao", "conv_tol": 1e-7},
+ conv_tol=1e-7,
+ gradient_kwargs={"eri_contraction": "direct", "conv_tol": 1e-7},
)
# The default penalty uses the smoothed Levine--Coe--Martinez form, the same
# formulation as geomeTRIC's built-in conical-intersection engine.
@@ -89,13 +91,14 @@ def energy_and_gradient(mol_at_step):
method = as_pyscf_method(mol, energy_and_gradient)
mol_ci = geometric_solver.optimize(
method,
- maxsteps=40,
+ maxsteps=100,
convergence_grms=1e-4,
convergence_gmax=2e-4,
)
- print_geometry("MECI geometry", mol_ci)
+ print_geometry("S0/S1 MECP geometry", mol_ci)
- # Report the final two surfaces at the located geometry.
+ # Report the final two surfaces (MP2 ground S0 + tracked excited S1) at the
+ # located geometry.
(e_lo, _), (e_hi, _) = paired(mol_ci.atom_coords(unit="Bohr"))
- print(f"\nFinal surface energies: lower = {e_lo:.10f} Eh, "
- f"upper = {e_hi:.10f} Eh, gap = {abs(e_hi - e_lo):.2e} Eh")
+ print(f"\nFinal surface energies: S0 = {e_lo:.10f} Eh, "
+ f"S1 = {e_hi:.10f} Eh, gap = {abs(e_hi - e_lo):.2e} Eh")
From 73e5f6cfe93206d2317c558441599bd8e7682d7f Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Mon, 22 Jun 2026 15:47:57 +0200
Subject: [PATCH 4/9] Switch MECP example to spin-flip ADC triplet MECI
Replace the ground/excited MECP setup (which hit SCF convergence trouble near
the diradical seam) with a spin-flip ADC(2) calculation from a triplet
unrestricted Hartree-Fock reference on twisted ethylene. The first two
spin-flip states recover the singlet ground state S0 and the first excited
singlet S1, so the S0/S1 crossing is optimised as a single-solve
excited/excited MECI -- no separate ground-state surface, no derivative
couplings. Verified the optimisation runs all steps cleanly (the MECP crashed).
- pyscf_adcc_mecp: triplet UHF reference (spin=2), n_spin_flip=4, states=(0,1).
- gradients.rst: note the spin-flip ADC route to S0/S1 conical intersections.
---
docs/gradients.rst | 12 ++++-
examples/optimization/pyscf_adcc_mecp.py | 56 +++++++++++++++---------
2 files changed, 46 insertions(+), 22 deletions(-)
diff --git a/docs/gradients.rst b/docs/gradients.rst
index 5dac98cb..068ff90f 100644
--- a/docs/gradients.rst
+++ b/docs/gradients.rst
@@ -146,8 +146,16 @@ smoothing::
method = as_pyscf_method(mol, lambda m: objective(m.atom_coords(unit="Bohr")))
mol_ci = geometric_solver.optimize(method, maxsteps=20)
-A complete worked example is provided in
-``examples/optimization/pyscf_adcc_mecp.py``.
+**Spin-flip ADC for S0/S1 conical intersections.** Because both surfaces of a
+MECI pair come from a single ADC solve, the cleanest route to an S0/S1 crossing
+is a *spin-flip* ADC calculation from a triplet unrestricted Hartree--Fock
+reference: the first spin-flip state recovers the closed-shell singlet ground
+state ``S0`` and the next one the first excited singlet ``S1``, so the S0/S1 CI
+is optimised as an excited/excited MECI (no separate ground-state surface, which
+avoids the convergence problems of a ground/excited MECP near a diradical
+seam). Use ``n_spin_flip`` in place of ``n_singlets`` for these scans. The
+worked example in ``examples/optimization/pyscf_adcc_mecp.py`` follows this
+spin-flip setup on twisted ethylene.
The two-electron term can be evaluated with two different strategies, selected
through the ``eri_contraction`` keyword (see
diff --git a/examples/optimization/pyscf_adcc_mecp.py b/examples/optimization/pyscf_adcc_mecp.py
index 25a732d9..99b98580 100644
--- a/examples/optimization/pyscf_adcc_mecp.py
+++ b/examples/optimization/pyscf_adcc_mecp.py
@@ -1,14 +1,23 @@
#!/usr/bin/env python3
-"""Minimum-energy crossing point (MECP) of the S0/S1 seam with adcc gradients.
-
-This example locates the minimum on the crossing seam between the ground state
-S0 and the first excited singlet S1 of twisted ethylene (C2H4) -- the textbook
-S0/S1 MECP / conical-intersection system -- using the paired-state scanner and
-the Levine--Coe--Martinez penalty objective. The
-:class:`adcc.PairedStateGradientScanner` evaluates the ground-state MP2 surface
-and one tracked ADC excited surface at a single geometry from a single SCF +
-single ADC, and :class:`adcc.MECPObjective` combines them into one
-``(energy, gradient)`` for geomeTRIC. **No derivative couplings** are required.
+"""Spin-flip ADC conical-intersection optimisation with adcc gradients.
+
+This example locates the S0/S1 conical intersection of twisted ethylene
+(C2H4), the textbook minimal MECI system, using a **spin-flip ADC** setup:
+
+* the Hartree--Fock reference is the **triplet** ground state (unrestricted
+ SCF, ``spin = 2``), the natural diradical reference for the twisted geometry;
+* a single spin-flip ADC(2) solve produces both the singlet ground state
+ ``S0`` (the first spin-flip state) and the first excited singlet ``S1``;
+* the :class:`adcc.PairedStateGradientScanner` evaluates these two surfaces at
+ one geometry from a single SCF + single ADC, and
+ :class:`adcc.MECPObjective` combines them into one ``(energy, gradient)``
+ for geomeTRIC's penalty conical-intersection driver.
+
+Because both surfaces come from the *same* spin-flip ADC solve, the
+S0/S1 crossing is treated as a proper excited/excited MECI -- no separate
+ground-state surface is optimised, which avoids the convergence problems of a
+ground/excited MECP formulation near the diradical region. **No derivative
+couplings** are required.
geomeTRIC is optional; install it separately, e.g. ``pip install geometric`` or
``pip install adcc[geomopt]`` once the optional extra is available.
@@ -21,8 +30,14 @@
from pyscf.geomopt import as_pyscf_method, geometric_solver
-def twisted_ethylene(twist_deg=80.0, basis="sto-3g"):
- """Build a twisted ethylene PySCF SCF near the S1/S0 crossing seam."""
+def twisted_ethylene(twist_deg=90.0, basis="6-31g"):
+ """Build an unrestricted triplet-HF ethylene near the S0/S1 crossing seam.
+
+ The triplet is the diradical reference for spin-flip ADC: the first
+ spin-flip state recovers the closed-shell singlet ground state ``S0`` and
+ the next one the first excited singlet ``S1``, so the S0/S1 crossing is
+ accessible as a single-solve excited/excited MECI.
+ """
tw = np.deg2rad(twist_deg)
c1 = np.array([0.0, 0.0, 0.0])
c2 = np.array([1.34, 0.0, 0.0])
@@ -46,9 +61,10 @@ def twist(p, angle, pivot):
)
mol = gto.M(
atom=atom_str, basis=basis, unit="Angstrom",
+ spin=2, # 2S = 2 -> triplet reference for spin-flip ADC
symmetry=False, verbose=0, parse_arg=False,
)
- mf = scf.RHF(mol)
+ mf = scf.UHF(mol)
mf.conv_tol = 1e-10
mf.conv_tol_grad = 1e-7
mf.max_cycles = 250
@@ -63,15 +79,16 @@ def print_geometry(title, mol):
if __name__ == "__main__":
- scfres = twisted_ethylene(twist_deg=80.0, basis="6-31g")
+ scfres = twisted_ethylene(twist_deg=90.0, basis="6-31g")
mol = scfres.mol
- print_geometry("Initial twisted geometry", mol)
+ print_geometry("Initial twisted geometry (triplet reference)", mol)
paired = adcc.PairedStateGradientScanner(
scfres,
method="adc2",
- lower="mp2", upper=0, # S0 (MP2 ground) + S1 (1st tracked excited)
- n_singlets=3,
+ states=(0, 1), # first two spin-flip states: S0 (singlet GS)
+ # and S1 (first excited singlet) -- a MECI
+ n_spin_flip=4, # spin-flip ADC, unrestricted reference only
follow="overlap",
conv_tol=1e-7,
gradient_kwargs={"eri_contraction": "direct", "conv_tol": 1e-7},
@@ -95,10 +112,9 @@ def energy_and_gradient(mol_at_step):
convergence_grms=1e-4,
convergence_gmax=2e-4,
)
- print_geometry("S0/S1 MECP geometry", mol_ci)
+ print_geometry("S0/S1 MECI geometry", mol_ci)
- # Report the final two surfaces (MP2 ground S0 + tracked excited S1) at the
- # located geometry.
+ # Report the final two spin-flip surfaces (S0 + S1) at the located geometry.
(e_lo, _), (e_hi, _) = paired(mol_ci.atom_coords(unit="Bohr"))
print(f"\nFinal surface energies: S0 = {e_lo:.10f} Eh, "
f"S1 = {e_hi:.10f} Eh, gap = {abs(e_hi - e_lo):.2e} Eh")
From 6b93dd9215cfb888fe0815d35916f730d0212385 Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Mon, 22 Jun 2026 16:33:13 +0200
Subject: [PATCH 5/9] Use follow="index" for MECI pair (density tracking fights
adiabatic seam)
The SF-ADC example diverged to a 90 mEh gap because follow="overlap" tried to
preserve fixed state character across geometries and, once S0/S1 densities
blur at the seam, flipped a slot onto a higher root. For a MECI the two surfaces
must be the adiabatically lowest roots at each geometry (the seam is defined by
an energy-ordered degeneracy), so follow="index" is correct -- exactly the
contract geomeTRIC's ConicalIntersection engine expects from its sub-engines.
- pyscf_adcc_mecp: follow="index", sigma=200, alpha=0.025, tighter geomeTRIC
grms/gmax; verified end-to-end convergence to gap = 5.9e-4 Eh (~0.016 eV).
- gradients.rst: note the follow="index" vs follow="overlap" distinction
for paired MECI optimisation.
---
docs/gradients.rst | 12 ++++++++++++
examples/optimization/pyscf_adcc_mecp.py | 25 ++++++++++++++++--------
2 files changed, 29 insertions(+), 8 deletions(-)
diff --git a/docs/gradients.rst b/docs/gradients.rst
index 068ff90f..f1bf8853 100644
--- a/docs/gradients.rst
+++ b/docs/gradients.rst
@@ -157,6 +157,18 @@ seam). Use ``n_spin_flip`` in place of ``n_singlets`` for these scans. The
worked example in ``examples/optimization/pyscf_adcc_mecp.py`` follows this
spin-flip setup on twisted ethylene.
+.. note::
+
+ For a MECI pair use ``follow="index"`` (the adiabatically lowest roots), not
+ ``follow="overlap"``. The crossing seam is defined by a degeneracy of the
+ *energy-ordered* pair, so the two tracked surfaces should simply be the
+ two lowest roots at each geometry; the single-surface density-overlap tracker
+ is meant to preserve a *fixed* state character across a geometry change --
+ exactly what a MECI optimiser should *not* do near the seam, where it would
+ fight the adiabatic reordering and push a slot onto a higher root. This
+ matches the contract geomeTRIC's ``ConicalIntersection`` engine expects from
+ its sub-engines ("roots 0 and 1, energy-sorted").
+
The two-electron term can be evaluated with two different strategies, selected
through the ``eri_contraction`` keyword (see
:ref:`gradients-eri-contraction` below). For the PySCF backend the
diff --git a/examples/optimization/pyscf_adcc_mecp.py b/examples/optimization/pyscf_adcc_mecp.py
index 99b98580..895062a3 100644
--- a/examples/optimization/pyscf_adcc_mecp.py
+++ b/examples/optimization/pyscf_adcc_mecp.py
@@ -89,13 +89,22 @@ def print_geometry(title, mol):
states=(0, 1), # first two spin-flip states: S0 (singlet GS)
# and S1 (first excited singlet) -- a MECI
n_spin_flip=4, # spin-flip ADC, unrestricted reference only
- follow="overlap",
- conv_tol=1e-7,
- gradient_kwargs={"eri_contraction": "direct", "conv_tol": 1e-7},
+ # For a MECI pair the two surfaces are the adiabatically lowest states
+ # at each geometry (the seam is defined by an energy-ordered degeneracy,
+ # not by fixed state character). Density-overlap tracking -- which
+ # serves single-surface optimization -- would *fight* the adiabatic
+ # reordering near the seam and flip a slot onto a higher root. Instead
+ # follow="index" returns positional roots (0, 1) and lets the scanner
+ # energy-sort them, exactly the contract geomeTRIC's conical-intersection
+ # engine expects from its sub-engines.
+ follow="index",
+ conv_tol=1e-8,
+ gradient_kwargs={"eri_contraction": "direct", "conv_tol": 1e-8},
)
# The default penalty uses the smoothed Levine--Coe--Martinez form, the same
- # formulation as geomeTRIC's built-in conical-intersection engine.
- objective = adcc.MECPObjective(paired)
+ # formulation as geomeTRIC's built-in conical-intersection engine. A larger
+ # sigma enforces the degeneracy harder for a tighter final gap.
+ objective = adcc.MECPObjective(paired, sigma=200.0, alpha=0.025)
def energy_and_gradient(mol_at_step):
energy, gradient = objective(mol_at_step.atom_coords(unit="Bohr"))
@@ -108,9 +117,9 @@ def energy_and_gradient(mol_at_step):
method = as_pyscf_method(mol, energy_and_gradient)
mol_ci = geometric_solver.optimize(
method,
- maxsteps=100,
- convergence_grms=1e-4,
- convergence_gmax=2e-4,
+ maxsteps=50,
+ convergence_grms=3e-5,
+ convergence_gmax=1e-4,
)
print_geometry("S0/S1 MECI geometry", mol_ci)
From 70f19e07742dcc22d9e147757959ad178bf3779e Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Mon, 22 Jun 2026 17:08:07 +0200
Subject: [PATCH 6/9] Warn when follow="overlap" is used for a MECI pair
Density-overlap tracking is designed for the single-surface scanner (follow one
fixed state character across geometries). For an excited/excited MECI pair the
seam is defined by a degeneracy of the energy-ordered two lowest roots, so
tracking by overlap fights the adiabatic reordering near the seam and can flip a
slot onto a higher root -- diverging the penalty optimisation (the failure mode
that left the SF-ADC example at a 90 mEh gap). Add a UserWarning at scanner
construction pointing at follow="index", which is the contract geomeTRIC's
ConicalIntersection engine expects from its sub-engines. The warning is
specific to the MECI pair form; ground/excited MECP only tracks one root, so
follow="overlap" stays a valid choice there.
- paired_scanner: warn on follow="overlap" + PairedExcitedStateTarget.
- geomopt_mecp_test: assert the warning fires for MECI, does NOT fire for MECP
or for follow="index"; silence the (deliberate) warning on the eight tests
that exercise the overlap-tracking mechanism itself.
---
adcc/gradients/paired_scanner.py | 26 +++++++++++++++++++
adcc/tests/geomopt_mecp_test.py | 44 ++++++++++++++++++++++++++++++++
2 files changed, 70 insertions(+)
diff --git a/adcc/gradients/paired_scanner.py b/adcc/gradients/paired_scanner.py
index 8c81413c..08153ab6 100644
--- a/adcc/gradients/paired_scanner.py
+++ b/adcc/gradients/paired_scanner.py
@@ -37,6 +37,7 @@
from __future__ import annotations
+import warnings
from dataclasses import dataclass, field
from typing import Any, Optional, Sequence
@@ -188,6 +189,31 @@ def __init__(self, scfres, *,
run_adc_kwargs,
)
+ # Density-overlap tracking is designed for the single-surface scanner
+ # (follow one fixed state character across geometries). For an
+ # excited/excited MECI pair the seam is defined by a degeneracy of the
+ # energy-ordered two lowest roots, not by fixed state character, so
+ # tracking by overlap fights the adiabatic reordering near the seam --
+ # once the two roots' densities blur together it can lock a slot onto a
+ # higher root and the optimisation diverges. Use follow="index" (the
+ # scanner energy-sorts the positional roots and the penalty drives the
+ # gap), which is the contract geomeTRIC's ConicalIntersection engine
+ # expects from its sub-engines.
+ if self.follow == "overlap" and isinstance(
+ self.paired_target, PairedExcitedStateTarget):
+ warnings.warn(
+ "follow='overlap' is not recommended for a MECI pair "
+ "(excited/excited): near the crossing seam the two tracked "
+ "roots' densities become near-degenerate and overlap tracking "
+ "can flip a slot onto a higher root, diverging the penalty "
+ "optimisation. Use follow='index' (the two lowest adiabatic "
+ "roots, energy-sorted) for MECI optimisations; this is also "
+ "what geomeTRIC's ConicalIntersection engine expects from its "
+ "sub-engines.",
+ UserWarning,
+ stacklevel=2,
+ )
+
# Plural tracking state. For MECP slot 0 (ground) is never tracked:
# its descriptor / index stay ``None`` forever.
self.previous_descriptors: tuple[Optional[_TrackingDescriptor],
diff --git a/adcc/tests/geomopt_mecp_test.py b/adcc/tests/geomopt_mecp_test.py
index ea9f1e72..778b9e0c 100644
--- a/adcc/tests/geomopt_mecp_test.py
+++ b/adcc/tests/geomopt_mecp_test.py
@@ -292,6 +292,7 @@ def test_paired_calc_new_returns_both_surfaces():
assert result["gradient_upper"].shape == (9,)
+@pytest.mark.filterwarnings("ignore::UserWarning")
def test_paired_excited_overlap_tracking_seeds_and_tracks_both_slots():
scanner = adcc.PairedStateGradientScanner(
_h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
@@ -358,6 +359,7 @@ def __init__(self, index, transition_dm, state_diffdm, omega):
self.excitation_energy = omega
+@pytest.mark.filterwarnings("ignore::UserWarning")
def test_distinctness_guard_keeps_slots_apart_under_collapse():
# Both previous descriptors would, independently, best-match candidate 0.
# The joint selection must keep the two roots distinct: one slot takes 0,
@@ -500,6 +502,42 @@ def test_paired_invalid_follow_raises_eagerly():
)
+def test_paired_meci_overlap_follow_emits_warning():
+ # follow="overlap" fights the adiabatic energy-ordering that defines a MECI
+ # seam and can flip a slot onto a higher root near degeneracy. The scanner
+ # warns (it does not raise -- the user may know what they are doing for
+ # well-separated states) and points at follow="index".
+ with pytest.warns(UserWarning, match="follow='overlap' is not recommended"):
+ adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="overlap",
+ )
+
+
+def test_paired_mecp_overlap_follow_does_not_warn():
+ # The warning is specific to the excited/excited MECI pair. A ground/excited
+ # MECP only tracks one (excited) root, so follow="overlap" is a legitimate
+ # choice there and must not warn.
+ import warnings as _w
+ with _w.catch_warnings():
+ _w.simplefilter("error", UserWarning)
+ adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", lower="mp2", upper=0, n_singlets=3,
+ follow="overlap",
+ )
+
+
+def test_paired_meci_index_follow_does_not_warn():
+ # The recommended setting for a MECI pair must not warn.
+ import warnings as _w
+ with _w.catch_warnings():
+ _w.simplefilter("error", UserWarning)
+ adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="index",
+ )
+
+
def test_paired_validates_coordinate_shape():
scanner = adcc.PairedStateGradientScanner(
_h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
@@ -665,6 +703,7 @@ def test_penalty_rejects_non_finite_gradients():
# Paired threshold enforcement (finding 3) for MECI and MECP overlap tracking.
# ---------------------------------------------------------------------------
+@pytest.mark.filterwarnings("ignore::UserWarning")
def test_paired_overlap_tracking_below_min_score_raises():
scanner = adcc.PairedStateGradientScanner(
_h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
@@ -676,6 +715,7 @@ def test_paired_overlap_tracking_below_min_score_raises():
scanner(scanner.initial_coords)
+@pytest.mark.filterwarnings("ignore::UserWarning")
def test_paired_overlap_tracking_ambiguous_gap_raises():
scanner = adcc.PairedStateGradientScanner(
_h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
@@ -713,6 +753,7 @@ def test_mecp_overlap_tracking_ambiguous_gap_raises():
# tracking_min_gap with only two states (finding 4).
# ---------------------------------------------------------------------------
+@pytest.mark.filterwarnings("ignore::UserWarning")
def test_tracking_min_gap_raises_when_only_two_excited_states():
# n_singlets == 2 leaves no per-channel alternative to compare a gap
# against; requesting ambiguity detection must surface that clearly.
@@ -726,6 +767,7 @@ def test_tracking_min_gap_raises_when_only_two_excited_states():
scanner(scanner.initial_coords)
+@pytest.mark.filterwarnings("ignore::UserWarning")
def test_tracking_min_gap_inert_at_two_states_with_default():
# With the default tracking_min_gap == 0 the two-state case does not raise;
# the per-channel gap is reported as inf (no measurable alternative).
@@ -771,6 +813,7 @@ def test_mecp_overlap_tracking_seeds_and_tracks_excited_root():
assert excited_tr.switched is False
+@pytest.mark.filterwarnings("ignore::UserWarning")
def test_paired_overlap_tracking_records_per_channel_diagnostics():
scanner = adcc.PairedStateGradientScanner(
_h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
@@ -947,6 +990,7 @@ def test_paired_meci_duplicate_index_mode_raises_distinctness():
# Joint tie-break toward the smaller excitation-energy gap (finding 13).
# ---------------------------------------------------------------------------
+@pytest.mark.filterwarnings("ignore::UserWarning")
def test_joint_selection_tie_break_prefers_smaller_excitation_gap(monkeypatch):
# Both slots track via overlap with two distinct ordered pairs tied on the
# summed overlap score but differing in excitation-energy gap; the joint
From 1a63b593f6d7ab546ca3d3ea9b03d72a69c3990b Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Mon, 22 Jun 2026 22:39:34 +0200
Subject: [PATCH 7/9] Keep MECP alpha=0 penalty gradient at the seam; switch
MECI smoke+docs to follow=index
---
adcc/gradients/mecp.py | 52 +++++++++++++------
adcc/tests/functionality_geomopt_mecp_test.py | 39 +++++++++++---
adcc/tests/geomopt_mecp_test.py | 43 +++++++++++----
docs/gradients.rst | 2 +-
4 files changed, 100 insertions(+), 36 deletions(-)
diff --git a/adcc/gradients/mecp.py b/adcc/gradients/mecp.py
index 9956e281..fc380d2b 100644
--- a/adcc/gradients/mecp.py
+++ b/adcc/gradients/mecp.py
@@ -86,13 +86,23 @@ def mecp_penalty(e_lower, g_lower, e_upper, g_upper, *, sigma=DEFAULT_SIGMA,
seam (``E_dif = 0``); the ``sigma`` weight controls how hard the degeneracy
is enforced versus minimising the average energy.
- At exact degeneracy (``e_upper == e_lower``) the penalty and its gradient
- vanish, leaving ``E_obj = E_avg`` and ``G_obj = G_avg`` -- the objective
- reduces to an average-surface optimisation pinned to the crossing. This
- includes the ``alpha == 0`` raw energy-difference mode: the limiting
- penalty ``sigma * E_dif`` tends to zero as ``E_dif -> 0`` even though the
- un-guarded expression ``sigma * E_dif**2 / (E_dif + 0)`` is ``0/0`` there,
- so the exact-degeneracy point is short-circuited explicitly below.
+ At exact degeneracy (``e_upper == e_lower``) the penalty energy vanishes.
+ For the smoothed mode (``alpha > 0``) the penalty *gradient* vanishes there
+ too, so the objective reduces to an average-surface optimisation
+ (``E_obj = E_avg``, ``G_obj = G_avg``) pinned to the crossing -- the LCM
+ objective is continuously extended to the seam.
+
+ The raw energy-difference mode (``alpha == 0``) is the singular limit of the
+ same formula and is evaluated in closed form as
+ ``E_pen = sigma * E_dif / n_states2`` and ``G_pen = sigma * G_dif /
+ n_states2`` -- with no division, so the sub-DBL_MIN underflow regime is
+ handled as well as the exact ``E_dif == 0`` point. In this mode the
+ penalty energy still vanishes at the seam, but its gradient tends to the
+ constant ``sigma * G_dif`` (the force that pins the optimiser to the
+ crossing) and is *not* dropped at degeneracy; this is the genuine
+ continuous extension, not a vanishing one. ``alpha == 0`` therefore gives
+ a harder, non-smooth penalty; ``alpha > 0`` smooths it to a vanishing
+ gradient exactly at the seam.
Parameters
----------
@@ -127,16 +137,24 @@ def mecp_penalty(e_lower, g_lower, e_upper, g_upper, *, sigma=DEFAULT_SIGMA,
g_avg = 0.5 * (g_lower + g_upper)
# Smoothed penalty (Levine--Coe--Martinez / geomeTRIC ConicalIntersection).
- # At exact degeneracy the penalty term and its gradient vanish by
- # construction (the LCM objective is continuously extended to the seam).
- # This is the only singular point: with alpha == 0 the raw
- # energy-difference mode would otherwise evaluate 0.0 / .0 here, so we
- # short-circuit the seam explicitly. The guard also covers the (only
- # theoretical) sub-DBL_MIN underflow where e_dif**2 flushes to zero while
- # alpha == 0.
- if e_dif == 0.0 or (alpha == 0.0 and abs(e_dif) < 1e-300):
- e_pen = 0.0
- g_pen_scale = 0.0
+ #
+ # For alpha > 0 the smoothed formula is well-defined everywhere (denom =
+ # E_dif + alpha >= alpha > 0); at exact degeneracy (E_dif == 0) the penalty
+ # energy *and its gradient* vanish by construction, so the objective
+ # reduces to the average surface -- the LCM objective is continuously
+ # extended to the seam.
+ #
+ # The raw energy-difference mode (alpha == 0) is the singular limit: the
+ # un-guarded expression sigma * E_dif**2 / E_dif would be 0/0 at the seam.
+ # Its genuine continuous extension is the linear penalty sigma * E_dif:
+ # the energy still vanishes at the seam, but the gradient tends to the
+ # constant sigma * G_dif -- the force that pins the optimiser to the
+ # crossing -- and must NOT be dropped at E_dif == 0. Evaluate this branch
+ # without any division so the sub-DBL_MIN underflow regime (E_dif**2
+ # flushing to zero in float64) is handled too, not just the exact point.
+ if alpha == 0.0:
+ e_pen = sigma * e_dif / _N_STATES2
+ g_pen_scale = sigma / _N_STATES2
else:
denom = e_dif + alpha
e_pen = sigma * (e_dif * e_dif) / (denom * _N_STATES2)
diff --git a/adcc/tests/functionality_geomopt_mecp_test.py b/adcc/tests/functionality_geomopt_mecp_test.py
index aa422c46..a7604422 100644
--- a/adcc/tests/functionality_geomopt_mecp_test.py
+++ b/adcc/tests/functionality_geomopt_mecp_test.py
@@ -94,9 +94,14 @@ def test_paired_scanner_objective_drives_meci_optimization():
from pyscf.geomopt import as_pyscf_method, geometric_solver
scfres = _ethylene_scf()
+ # follow="index" for a MECI pair (the recommended setting, matching the
+ # example and the docs note): the two lowest adiabatic roots, energy-sorted,
+ # the contract geomeTRIC's ConicalIntersection engine expects from its
+ # sub-engines. (follow="overlap" for a MECI pair is a documented footgun --
+ # see the PairedStateGradientScanner construction warning.)
scanner = adcc.PairedStateGradientScanner(
scfres, method="adc2", states=(0, 1), n_singlets=3,
- follow="overlap", conv_tol=1e-8,
+ follow="index", conv_tol=1e-8,
gradient_kwargs={"eri_contraction": "full_ao"},
)
objective = adcc.MECPObjective(scanner) # default LCM penalty
@@ -105,6 +110,11 @@ def test_paired_scanner_objective_drives_meci_optimization():
def energy_and_gradient(mol_at_step):
energy, gradient = objective(mol_at_step.atom_coords(unit="Bohr"))
+ # The penalty must hand geomeTRIC a finite objective every step; a
+ # diverging run (the regression that motivated the MECI follow='overlap'
+ # warning) surfaces here rather than as a silent blow-up.
+ assert np.isfinite(energy)
+ assert np.all(np.isfinite(gradient))
e_lo, e_hi = objective.last_pair[0][0], objective.last_pair[1][0]
seen_gaps.append(abs(e_hi - e_lo))
return energy, gradient
@@ -118,11 +128,15 @@ def energy_and_gradient(mol_at_step):
# Multiple steps ran and both surfaces were evaluated every step.
assert len(seen_gaps) >= 2
assert all(np.isfinite(gap) for gap in seen_gaps)
- assert scanner.last_trackings is not None
- # Tracking across steps kept the two followed roots distinct throughout.
- slot0, slot1 = scanner.last_trackings
- assert slot0 is not None and slot1 is not None
- assert slot0.index != slot1.index
+ # The penalty optimisation should drive the gap down (or at least not let
+ # it blow up) -- "driving toward the seam", not merely "did not crash". A
+ # generous factor keeps the assertion robust over a tiny 3-step line search
+ # while still catching a diverging run that doubles or triples the gap.
+ assert seen_gaps[-1] <= 2.0 * seen_gaps[0] + 1e-9
+ # follow="index" keeps the two lowest adiabatic roots structurally distinct
+ # (index mode records no per-step tracking diagnostic, by design); the
+ # energy-sorted pair stays apart.
+ assert objective.last_pair[0][0] <= objective.last_pair[1][0]
def test_mecp_objective_calc_new_contract_end_to_end():
@@ -155,8 +169,15 @@ def test_mecp_ground_excited_pair_drives_optimization():
)
objective = adcc.MECPObjective(scanner)
+ seen_gaps = []
+
def energy_and_gradient(mol_at_step):
- return objective(mol_at_step.atom_coords(unit="Bohr"))
+ energy, gradient = objective(mol_at_step.atom_coords(unit="Bohr"))
+ assert np.isfinite(energy)
+ assert np.all(np.isfinite(gradient))
+ e_lo, e_hi = objective.last_pair[0][0], objective.last_pair[1][0]
+ seen_gaps.append(abs(e_hi - e_lo))
+ return energy, gradient
method = as_pyscf_method(scfres.mol, energy_and_gradient)
# A couple of steps confirm the MECP objective drives geomeTRIC without
@@ -169,3 +190,7 @@ def energy_and_gradient(mol_at_step):
e_lo, e_hi = objective.last_pair[0][0], objective.last_pair[1][0]
assert np.isfinite(e_lo) and np.isfinite(e_hi)
assert e_lo <= e_hi
+ # Multiple steps ran and the gap did not blow up (driving toward the seam).
+ assert len(seen_gaps) >= 2
+ assert all(np.isfinite(gap) for gap in seen_gaps)
+ assert seen_gaps[-1] <= 2.0 * seen_gaps[0] + 1e-9
diff --git a/adcc/tests/geomopt_mecp_test.py b/adcc/tests/geomopt_mecp_test.py
index 778b9e0c..321f53a2 100644
--- a/adcc/tests/geomopt_mecp_test.py
+++ b/adcc/tests/geomopt_mecp_test.py
@@ -656,35 +656,56 @@ def test_penalty_matches_geometric_oracle_at_degeneracy():
# Exact-degeneracy guard at alpha == 0 (findings 1 & 5) and finiteness guards.
# ---------------------------------------------------------------------------
-def test_penalty_alpha_zero_exact_degeneracy_collapses_to_average():
- # At (alpha == 0, e_dif == 0) the raw energy-difference mode would
- # otherwise evaluate 0.0 / 0.0; the guard short-circuits to the average
- # surface so the penalty vanishes exactly, matching the docstring.
+def test_penalty_alpha_zero_exact_degeneracy_keeps_penalty_gradient():
+ # At (alpha == 0, e_dif == 0) the raw energy-difference mode is evaluated in
+ # closed form (no division): the penalty energy vanishes, but its gradient
+ # tends to the constant sigma * G_dif -- the force that pins the optimiser
+ # to the crossing -- which must NOT be dropped at the seam. The objective
+ # gradient is therefore G_avg + sigma * G_dif, the genuine continuous
+ # extension, not the bare average surface.
g_lo = np.array([[1.0, 0.0, 0.0]])
g_hi = np.array([[0.0, 2.0, 0.0]])
e = -1.2
- energy, gradient = mecp_penalty(e, g_lo, e, g_hi, sigma=2.0, alpha=0.0)
- assert energy == pytest.approx(e)
- assert_allclose(gradient, 0.5 * (g_lo + g_hi))
+ sigma = 2.0
+ energy, gradient = mecp_penalty(e, g_lo, e, g_hi, sigma=sigma, alpha=0.0)
+ assert energy == pytest.approx(e) # penalty energy vanishes at the seam
+ assert_allclose(gradient, 0.5 * (g_lo + g_hi) + sigma * (g_hi - g_lo))
def test_penalty_alpha_zero_near_degeneracy_is_finite():
- # The sub-DBL_MIN underflow guard (-- abs(e_dif) < 1e-300 with alpha == 0)
- # must also short-circuit rather than produce 0.0 / 0.0.
- tiny = 1e-310 # e_dif ** 2 flushes to zero in float64
+ # The alpha == 0 raw mode carries no division, so the whole sub-DBL_MIN
+ # underflow regime (E_dif**2 flushing to zero in float64) -- not just the
+ # exact E_dif == 0 point -- is well-defined rather than 0.0 / 0.0.
+ tiny = 1e-310 # E_dif ** 2 flushes to zero in float64
energy, gradient = mecp_penalty(
-1.2, np.zeros(3), -1.2 + tiny, np.zeros(3), sigma=2.0, alpha=0.0,
)
assert np.isfinite(energy)
+ assert_allclose(gradient, np.zeros(3))
def efun(off):
return mecp_penalty(-1.2, np.zeros(3), -1.2 + off, np.zeros(3),
sigma=2.0, alpha=0.0)[0]
- # The limit from outside the guard (alpha == 0) is the average as off -> 0.
+ # The limit from outside the exact point (alpha == 0) is the average as
+ # off -> 0 (the penalty energy sigma * off vanishes).
assert efun(tiny) == pytest.approx(-1.2, abs=1e-6)
+def test_penalty_alpha_zero_underflow_regime_is_finite():
+ # Regression for the former ZeroDivisionError gap: with the old short-
+ # circuit (abs(E_dif) < 1e-300) a value of E_dif in [1e-300, ~7e-162) still
+ # made denom ** 2 underflow to zero while alpha == 0, raising 0.0 / 0.0.
+ # The no-division raw mode handles the whole regime. 1e-200 sits squarely
+ # in the former gap (E_dif ** 2 flushes to zero, but no division happens).
+ tiny = 1e-200
+ energy, gradient = mecp_penalty(
+ -1.2, np.zeros(3), -1.2 + tiny, np.zeros(3), sigma=2.0, alpha=0.0,
+ )
+ assert np.isfinite(energy)
+ assert_allclose(gradient, np.zeros(3))
+
+
def test_penalty_rejects_non_finite_energies():
g = np.zeros(3)
with pytest.raises(RuntimeError, match="non-finite"):
diff --git a/docs/gradients.rst b/docs/gradients.rst
index f1bf8853..6f9517e5 100644
--- a/docs/gradients.rst
+++ b/docs/gradients.rst
@@ -117,7 +117,7 @@ convenience::
method="adc2",
states=(0, 1), # two tracked excited states (MECI)
n_singlets=5,
- follow="overlap",
+ follow="index", # two lowest adiabatic roots, energy-sorted
)
(e_lower, g_lower), (e_upper, g_upper) = paired(mol.atom_coords())
From e151640803a6adb50dc82a84da3f8a9d601ab8c8 Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Tue, 23 Jun 2026 15:59:24 +0200
Subject: [PATCH 8/9] Let scanners/objective accept a PySCF Mole directly; add
step_callback
---
adcc/gradients/mecp.py | 9 +++
adcc/gradients/scanner.py | 19 ++++++-
adcc/tests/geomopt_mecp_test.py | 57 +++++++++++++++++++
adcc/tests/geomopt_scanner_test.py | 25 ++++++++
docs/gradients.rst | 9 ++-
examples/optimization/pyscf_adcc_geoopt.py | 9 +--
.../optimization/pyscf_adcc_geoopt_excited.py | 21 ++++---
examples/optimization/pyscf_adcc_mecp.py | 11 ++--
8 files changed, 140 insertions(+), 20 deletions(-)
diff --git a/adcc/gradients/mecp.py b/adcc/gradients/mecp.py
index fc380d2b..24ac3e0d 100644
--- a/adcc/gradients/mecp.py
+++ b/adcc/gradients/mecp.py
@@ -217,6 +217,13 @@ def __init__(self, scanner, *, sigma=DEFAULT_SIGMA, alpha=DEFAULT_ALPHA,
self.last_gradient = None
self.last_pair = None # ((e_lower, g_lower), (e_upper, g_upper))
+ # Optional opt-in per-step hook, invoked as ``cb(energy, gradient)`` at
+ # the end of every ``__call__``. Lets examples/tests attach a one-line
+ # "is it converging?" printer (penalty energy / gradient norm / seam gap,
+ # which is readable off ``self.last_pair``) without a hand-rolled
+ # wrapping function around the objective.
+ self.step_callback: Optional[Callable] = None
+
def __call__(self, coords):
"""Return ``(energy, gradient)`` for Cartesian coordinates in Bohr."""
(e_lower, g_lower), (e_upper, g_upper) = self.scanner(coords)
@@ -227,6 +234,8 @@ def __call__(self, coords):
)
self.last_energy = e_obj
self.last_gradient = np.asarray(g_obj)
+ if self.step_callback is not None:
+ self.step_callback(e_obj, self.last_gradient)
return e_obj, np.asarray(g_obj)
def calc_new(self, coords):
diff --git a/adcc/gradients/scanner.py b/adcc/gradients/scanner.py
index 451853d7..cbdc6a11 100644
--- a/adcc/gradients/scanner.py
+++ b/adcc/gradients/scanner.py
@@ -31,7 +31,7 @@
from __future__ import annotations
from dataclasses import dataclass, field
-from typing import Any, Optional
+from typing import Any, Callable, Optional
import warnings
import numpy as np
@@ -137,6 +137,12 @@ def __init__(self, scfres, *,
self.last_gradient = None
self.last_tracking: Optional[TrackingResult] = None
+ # Optional opt-in per-step hook, invoked as ``cb(energy, gradient)`` at
+ # the end of every ``__call__``. Lets examples/tests attach a one-line
+ # "is it converging?" printer (energy / gradient norm / seam gap) without
+ # a hand-rolled wrapping function around the scanner.
+ self.step_callback: Optional[Callable] = None
+
@property
def natoms(self) -> int:
return len(self.atom_symbols)
@@ -163,6 +169,9 @@ def __call__(self, coords):
self.previous_index = self.last_tracking.index
else:
self.previous_index = self.target.state_index
+
+ if self.step_callback is not None:
+ self.step_callback(float(energy), np.asarray(grad.total))
return float(energy), np.asarray(grad.total)
def calc_new(self, coords):
@@ -198,6 +207,14 @@ def _normalise_target(self, target, method, state_index, mp_level,
)
def _coords_array(self, coords):
+ # Accept a PySCF ``Mole`` (or anything exposing ``atom_coords``) so the
+ # scanner can be plugged straight into ``as_pyscf_method`` /
+ # ``geometric_solver.optimize`` without a hand-rolled wrapper extracting
+ # Bohr coordinates. The array path below is unchanged for callers that
+ # already pass a Cartesian array.
+ ac = getattr(coords, "atom_coords", None)
+ if callable(ac) and not isinstance(coords, np.ndarray):
+ coords = ac(unit="Bohr")
coords = np.asarray(coords, dtype=float)
if coords.shape == (3 * self.natoms,):
coords = coords.reshape(self.natoms, 3)
diff --git a/adcc/tests/geomopt_mecp_test.py b/adcc/tests/geomopt_mecp_test.py
index 321f53a2..e924fa48 100644
--- a/adcc/tests/geomopt_mecp_test.py
+++ b/adcc/tests/geomopt_mecp_test.py
@@ -238,6 +238,63 @@ def __call__(self, coords):
# Paired-engine mechanics (one SCF + one ADC, two surfaces).
# ---------------------------------------------------------------------------
+def test_scanner_accepts_pyscf_mole_directly():
+ # The scanner reads ``atom_coords(unit="Bohr")`` on a PySCF ``Mole``, so it
+ # plugs straight into ``as_pyscf_method`` without a hand-rolled wrapper.
+ scfres = _h2o_scf()
+ scanner = adcc.PairedStateGradientScanner(
+ scfres, method="adc2", states=(0, 1), n_singlets=3,
+ follow="index", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ pair_via_mole = scanner(scfres.mol)
+ pair_via_coords = scanner(scanner.initial_coords)
+ # Passing the Mole must give the same result as passing its Bohr coords.
+ assert pair_via_mole[0][0] == pytest.approx(pair_via_coords[0][0], abs=1e-10)
+ assert pair_via_mole[1][0] == pytest.approx(pair_via_coords[1][0], abs=1e-10)
+ assert_allclose(pair_via_mole[0][1], pair_via_coords[0][1], atol=1e-10)
+
+
+def test_mecp_objective_accepts_pyscf_mole_directly():
+ # The objective forwards the Mole to its paired scanner, so it too can be
+ # passed straight to ``as_pyscf_method``.
+ scfres = _h2o_scf()
+ scanner = adcc.PairedStateGradientScanner(
+ scfres, method="adc2", states=(0, 1), n_singlets=3,
+ follow="index", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ objective = adcc.MECPObjective(scanner)
+ e_via_mole, g_via_mole = objective(scfres.mol)
+ e_via_coords, g_via_coords = objective(scanner.initial_coords)
+ assert e_via_mole == pytest.approx(e_via_coords, abs=1e-10)
+ assert_allclose(g_via_mole, g_via_coords, atol=1e-10)
+
+
+def test_mecp_objective_step_callback_is_invoked():
+ # The opt-in per-step hook receives the objective's (energy, gradient) and
+ # runs *after* last_pair is populated, so it can read the seam gap.
+ scanner = adcc.PairedStateGradientScanner(
+ _h2o_scf(), method="adc2", states=(0, 1), n_singlets=3,
+ follow="index", conv_tol=1e-9,
+ gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ objective = adcc.MECPObjective(scanner)
+ seen = []
+
+ def _cb(e, g):
+ # last_pair is already set when the callback fires.
+ lo, hi = objective.last_pair
+ seen.append((e, float(np.asarray(g).ravel()[0]), abs(hi[0] - lo[0])))
+
+ objective.step_callback = _cb
+ objective(scanner.initial_coords)
+ assert len(seen) == 1
+ assert seen[0][0] == pytest.approx(objective.last_energy, abs=1e-12)
+ assert np.isfinite(seen[0][1])
+ assert np.isfinite(seen[0][2])
+
+
def _assert_independent_match(scanner):
states = adcc.run_adc(scanner.last_scf,
**scanner.paired_target.kwargs())
diff --git a/adcc/tests/geomopt_scanner_test.py b/adcc/tests/geomopt_scanner_test.py
index 76603e62..0806b697 100644
--- a/adcc/tests/geomopt_scanner_test.py
+++ b/adcc/tests/geomopt_scanner_test.py
@@ -102,6 +102,31 @@ def test_calc_new_returns_geometric_engine_shape():
assert result["gradient"].shape == (9,)
+def test_scanner_accepts_pyscf_mole_directly():
+ # The scanner reads atom_coords(unit="Bohr") on a PySCF Mole, so it can be
+ # passed straight to as_pyscf_method without a hand-rolled wrapper.
+ scfres = _h2o_scf()
+ scanner = adcc.NuclearGradientScanner(
+ scfres, method="mp2", gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ e_via_mole, g_via_mole = scanner(scfres.mol)
+ e_via_coords, g_via_coords = scanner(scanner.initial_coords)
+ assert e_via_mole == pytest.approx(e_via_coords, abs=1e-10)
+ assert_allclose(g_via_mole, g_via_coords, atol=1e-10)
+
+
+def test_scanner_step_callback_is_invoked_with_energy_and_gradient():
+ scanner = adcc.NuclearGradientScanner(
+ _h2o_scf(), method="mp2", gradient_kwargs={"eri_contraction": "full_ao"},
+ )
+ seen = []
+ scanner.step_callback = lambda e, g: seen.append((e, float(np.asarray(g).sum())))
+ e, g = scanner(scanner.initial_coords)
+ assert len(seen) == 1
+ assert seen[0][0] == pytest.approx(e, abs=1e-12)
+ assert seen[0][1] == pytest.approx(float(np.asarray(g).sum()), abs=1e-10)
+
+
def test_run_adc_kwargs_are_forwarded_with_native_names():
scanner = adcc.NuclearGradientScanner(
_h2o_scf(), method="adc2", n_singlets=2, conv_tol=1e-7,
diff --git a/docs/gradients.rst b/docs/gradients.rst
index 6f9517e5..69780f7b 100644
--- a/docs/gradients.rst
+++ b/docs/gradients.rst
@@ -143,9 +143,16 @@ smoothing::
from pyscf.geomopt import as_pyscf_method, geometric_solver
objective = adcc.MECPObjective(paired)
- method = as_pyscf_method(mol, lambda m: objective(m.atom_coords(unit="Bohr")))
+ method = as_pyscf_method(mol, objective)
mol_ci = geometric_solver.optimize(method, maxsteps=20)
+The scanner and objective accept a PySCF ``Mole`` directly (they read
+``atom_coords(unit="Bohr")``), so they plug straight into
+:func:`pyscf.geomopt.as_pyscf_method` without a hand-rolled wrapper. For an
+optional per-step "is it converging?" printout set a ``step_callback`` hook
+(``cb(energy, gradient)``); the worked examples in
+``examples/optimization/`` show both patterns.
+
**Spin-flip ADC for S0/S1 conical intersections.** Because both surfaces of a
MECI pair come from a single ADC solve, the cleanest route to an S0/S1 crossing
is a *spin-flip* ADC calculation from a triplet unrestricted Hartree--Fock
diff --git a/examples/optimization/pyscf_adcc_geoopt.py b/examples/optimization/pyscf_adcc_geoopt.py
index 7a9a477c..92506136 100644
--- a/examples/optimization/pyscf_adcc_geoopt.py
+++ b/examples/optimization/pyscf_adcc_geoopt.py
@@ -29,13 +29,10 @@
scanner = adcc.NuclearGradientScanner(mf, method="mp2")
-def energy_and_gradient(mol_at_step):
- """PySCF geomopt callback: return energy and gradient for a Mole."""
- return scanner(mol_at_step.atom_coords(unit="Bohr"))
-
-
if __name__ == "__main__":
- method = as_pyscf_method(mol, energy_and_gradient)
+ # The scanner accepts a PySCF Mole directly (it reads atom_coords in Bohr),
+ # so it can be passed straight to as_pyscf_method without a wrapper.
+ method = as_pyscf_method(mol, scanner)
mol_eq = geometric_solver.optimize(method, maxsteps=20)
print("\nOptimised geometry (Bohr):")
diff --git a/examples/optimization/pyscf_adcc_geoopt_excited.py b/examples/optimization/pyscf_adcc_geoopt_excited.py
index d187f96e..22266ea6 100644
--- a/examples/optimization/pyscf_adcc_geoopt_excited.py
+++ b/examples/optimization/pyscf_adcc_geoopt_excited.py
@@ -33,10 +33,14 @@ def print_geometry(title, mol):
f"{xyz[0]:16.10f} {xyz[1]:16.10f} {xyz[2]:16.10f}")
-def make_energy_gradient(scanner):
- """PySCF geomopt callback: return energy and gradient for a Mole."""
- def energy_and_gradient(mol_at_step):
- energy, gradient = scanner(mol_at_step.atom_coords(unit="Bohr"))
+def report_step(scanner):
+ """Optional per-step "is it converging?" printer attached as a callback.
+
+ Implemented as a closure over the scanner so it can read the just-evaluated
+ excitation (``scanner.last_excitation``); the scanner invokes it as
+ ``cb(energy, gradient)`` at the end of every ``__call__``.
+ """
+ def _cb(energy, gradient):
excitation = scanner.last_excitation
state_info = ""
if excitation is not None:
@@ -49,8 +53,7 @@ def energy_and_gradient(mol_at_step):
f"{float((gradient ** 2).sum() ** 0.5):.6e} Eh/Bohr"
f"{state_info}"
)
- return energy, gradient
- return energy_and_gradient
+ return _cb
# Use Bohr coordinates and keep molecular symmetry disabled so the optimizer's
@@ -87,7 +90,8 @@ def energy_and_gradient(mol_at_step):
method="mp2",
gradient_kwargs={"conv_tol": 1e-7},
)
- mp2_method = as_pyscf_method(mol, make_energy_gradient(mp2_scanner))
+ mp2_scanner.step_callback = report_step(mp2_scanner)
+ mp2_method = as_pyscf_method(mol, mp2_scanner)
mol_mp2 = geometric_solver.optimize(
mp2_method,
maxsteps=20,
@@ -108,7 +112,8 @@ def energy_and_gradient(mol_at_step):
conv_tol=1e-6,
gradient_kwargs={"conv_tol": 1e-7},
)
- adc_method = as_pyscf_method(mol_mp2, make_energy_gradient(adc_scanner))
+ adc_scanner.step_callback = report_step(adc_scanner)
+ adc_method = as_pyscf_method(mol_mp2, adc_scanner)
mol_adc = geometric_solver.optimize(
adc_method,
maxsteps=20,
diff --git a/examples/optimization/pyscf_adcc_mecp.py b/examples/optimization/pyscf_adcc_mecp.py
index 895062a3..a03ad614 100644
--- a/examples/optimization/pyscf_adcc_mecp.py
+++ b/examples/optimization/pyscf_adcc_mecp.py
@@ -106,15 +106,18 @@ def print_geometry(title, mol):
# sigma enforces the degeneracy harder for a tighter final gap.
objective = adcc.MECPObjective(paired, sigma=200.0, alpha=0.025)
- def energy_and_gradient(mol_at_step):
- energy, gradient = objective(mol_at_step.atom_coords(unit="Bohr"))
+ # Optional per-step "is it converging?" printer attached as a callback; read
+ # the seam gap off objective.last_pair. The objective accepts a PySCF Mole
+ # directly (forwarded to the paired scanner), so it plugs straight into
+ # as_pyscf_method without a wrapper.
+ def _print_step(energy, gradient):
e_lo, e_hi = objective.last_pair[0][0], objective.last_pair[1][0]
gnorm = float((gradient ** 2).sum() ** 0.5)
print(f"E_pen = {energy:.12f} Eh, |g| = {gnorm:.6e} Eh/Bohr, "
f"gap = {abs(e_hi - e_lo):.6e} Eh")
- return energy, gradient
+ objective.step_callback = _print_step
- method = as_pyscf_method(mol, energy_and_gradient)
+ method = as_pyscf_method(mol, objective)
mol_ci = geometric_solver.optimize(
method,
maxsteps=50,
From 85a0de10990267bf2a61dbd9157d04004f2b3857 Mon Sep 17 00:00:00 2001
From: Maximilian Scheurer <11406414+maxscheurer@users.noreply.github.com>
Date: Wed, 24 Jun 2026 16:24:08 +0200
Subject: [PATCH 9/9] appease flake8
---
adcc/gradients/__init__.py | 110 +++++++++++++--------
adcc/tests/geomopt_scanner_test.py | 151 +++++++++++++++++++++--------
2 files changed, 176 insertions(+), 85 deletions(-)
diff --git a/adcc/gradients/__init__.py b/adcc/gradients/__init__.py
index 870ba5f6..f80601c6 100644
--- a/adcc/gradients/__init__.py
+++ b/adcc/gradients/__init__.py
@@ -33,21 +33,35 @@
from adcc.NParticleOperator import OperatorSymmetry, product_trace
from .TwoParticleDensityMatrix import TwoParticleDensityMatrix
from .orbital_response import (
- orbital_response, orbital_response_rhs, energy_weighted_density_matrix
+ orbital_response,
+ orbital_response_rhs,
+ energy_weighted_density_matrix,
)
from .amplitude_response import amplitude_relaxed_densities
-from .scanner import (ExcitedStateTarget, GroundStateTarget,
- NuclearGradientScanner, density_overlap_score)
-from .paired_scanner import (PairedExcitedStateTarget,
- PairedGroundExcitedStateTarget,
- PairedStateGradientScanner)
+from .scanner import (
+ ExcitedStateTarget,
+ GroundStateTarget,
+ NuclearGradientScanner,
+ density_overlap_score,
+)
+from .paired_scanner import (
+ PairedExcitedStateTarget,
+ PairedGroundExcitedStateTarget,
+ PairedStateGradientScanner,
+)
from .mecp import mecp_penalty, MECPObjective
__all__ = [
- "nuclear_gradient", "NuclearGradientScanner",
- "GroundStateTarget", "ExcitedStateTarget", "density_overlap_score",
- "PairedStateGradientScanner", "PairedExcitedStateTarget",
- "PairedGroundExcitedStateTarget", "mecp_penalty", "MECPObjective",
+ "nuclear_gradient",
+ "NuclearGradientScanner",
+ "GroundStateTarget",
+ "ExcitedStateTarget",
+ "density_overlap_score",
+ "PairedStateGradientScanner",
+ "PairedExcitedStateTarget",
+ "PairedGroundExcitedStateTarget",
+ "mecp_penalty",
+ "MECPObjective",
]
@@ -95,11 +109,9 @@ def _energy(self):
"""Compute energy based on density matrices
for testing purposes"""
if self.g1a is None:
- raise ValueError("No unrelaxed one-particle "
- "density available.")
+ raise ValueError("No unrelaxed one-particle density available.")
if self.g2a is None:
- raise ValueError("No unrelaxed two-particle "
- "density available.")
+ raise ValueError("No unrelaxed two-particle density available.")
ret = 0.0
hf = self.reference_state
for b in self.g1a.blocks_nonzero:
@@ -117,8 +129,7 @@ def dipole_moment_relaxed(self):
def dipole_moment_unrelaxed(self):
"""Returns the unrelaxed electric dipole moment"""
if self.g1a is None:
- raise ValueError("No unrelaxed one-particle "
- "density available.")
+ raise ValueError("No unrelaxed one-particle density available.")
hf = self.reference_state
return self.__dipole_moment_electric(self.g1a + hf.density)
@@ -129,22 +140,26 @@ def total(self):
def __dipole_moment_electric(self, dm):
dips = self.reference_state.operators.electric_dipole
- elec_dip = -1.0 * np.array(
- [product_trace(dm, dip) for dip in dips]
- )
+ elec_dip = -1.0 * np.array([product_trace(dm, dip) for dip in dips])
return elec_dip + self.reference_state.nuclear_dipole
-def nuclear_gradient(excitation_or_mp, conv_tol=1e-9, eri_contraction=None,
- eri_shell_chunk_size=1, eri_pair_chunk_size=None,
- eri_pair_density_storage="memory"):
+def nuclear_gradient(
+ excitation_or_mp,
+ conv_tol=1e-9,
+ eri_contraction=None,
+ eri_shell_chunk_size=1,
+ eri_pair_chunk_size=None,
+ eri_pair_density_storage="memory",
+):
if isinstance(excitation_or_mp, LazyMp):
mp = excitation_or_mp
elif isinstance(excitation_or_mp, Excitation):
mp = excitation_or_mp.ground_state
else:
- raise TypeError("Gradient can only be computed for "
- "Excitation or LazyMp object.")
+ raise TypeError(
+ "Gradient can only be computed for Excitation or LazyMp object."
+ )
timer = Timer()
hf = mp.reference_state
@@ -163,6 +178,7 @@ def nuclear_gradient(excitation_or_mp, conv_tol=1e-9, eri_contraction=None,
# For NOSYMMETRY operators (CVS-ADC1+), explicitly set transpose blocks
if g1o.symmetry == OperatorSymmetry.NOSYMMETRY:
from adcc.functions import transpose
+
g1o.vo = transpose(g1o.ov)
g1o.vc = transpose(g1o.cv)
if not g1o.is_zero_block("o2o1"):
@@ -183,16 +199,17 @@ def nuclear_gradient(excitation_or_mp, conv_tol=1e-9, eri_contraction=None,
if hf.has_core_occupied_space:
delta_IJ = hf.density.cc
- g2_hf.oooo = 0.25 * (- einsum("li,jk->ijkl", delta_ij, delta_ij)
- + einsum("ki,jl->ijkl", delta_ij, delta_ij))
+ g2_hf.oooo = 0.25 * (
+ -einsum("li,jk->ijkl", delta_ij, delta_ij)
+ + einsum("ki,jl->ijkl", delta_ij, delta_ij)
+ )
g2_hf.cccc = -0.5 * einsum("IK,JL->IJKL", delta_IJ, delta_IJ)
g2_hf.ococ = -1.0 * einsum("ik,JL->iJkL", delta_ij, delta_IJ)
g2_oresp.cccc = einsum("IK,JL->IJKL", delta_IJ, g1o.cc + delta_IJ)
- g2_oresp.ococ = (
- + einsum("ik,JL->iJkL", delta_ij, g1o.cc + 2.0 * delta_IJ)
- + einsum("ik,JL->iJkL", g1o.oo, delta_IJ)
- )
+ g2_oresp.ococ = +einsum(
+ "ik,JL->iJkL", delta_ij, g1o.cc + 2.0 * delta_IJ
+ ) + einsum("ik,JL->iJkL", g1o.oo, delta_IJ)
g2_oresp.oooo = einsum("ij,kl->kilj", delta_ij, g1o.oo)
g2_oresp.ovov = einsum("ij,ab->iajb", delta_ij, g1o.vv)
g2_oresp.cvcv = einsum("IJ,ab->IaJb", delta_IJ, g1o.vv)
@@ -211,13 +228,16 @@ def nuclear_gradient(excitation_or_mp, conv_tol=1e-9, eri_contraction=None,
g2_total = evaluate(g2_hf + g2a + g2_oresp)
else:
- g2_hf.oooo = 0.25 * (- einsum("li,jk->ijkl", delta_ij, delta_ij)
- + einsum("ki,jl->ijkl", delta_ij, delta_ij))
+ g2_hf.oooo = 0.25 * (
+ -einsum("li,jk->ijkl", delta_ij, delta_ij)
+ + einsum("ki,jl->ijkl", delta_ij, delta_ij)
+ )
g2_oresp.oooo = einsum("ij,kl->kilj", delta_ij, g1o.oo)
g2_oresp.ovov = einsum("ij,ab->iajb", delta_ij, g1o.vv)
- g2_oresp.ooov = (- einsum("kj,ia->ijka", delta_ij, g1o.ov)
- + einsum("ki,ja->ijka", delta_ij, g1o.ov))
+ g2_oresp.ooov = -einsum("kj,ia->ijka", delta_ij, g1o.ov) + einsum(
+ "ki,ja->ijka", delta_ij, g1o.ov
+ )
# scale for contraction with integrals
g2a.oovv *= 0.5
@@ -227,8 +247,9 @@ def nuclear_gradient(excitation_or_mp, conv_tol=1e-9, eri_contraction=None,
provider = hf.gradient_provider
if eri_contraction is None:
- if getattr(provider, "backend", None) == "pyscf" \
- and hasattr(provider, "correlated_gradient_direct"):
+ if getattr(provider, "backend", None) == "pyscf" and hasattr(
+ provider, "correlated_gradient_direct"
+ ):
eri_contraction = "direct"
else:
eri_contraction = "full_ao"
@@ -239,7 +260,8 @@ def nuclear_gradient(excitation_or_mp, conv_tol=1e-9, eri_contraction=None,
f"{sorted(valid_eri_contractions)}."
)
if eri_contraction == "direct" and not hasattr(
- provider, "correlated_gradient_direct"):
+ provider, "correlated_gradient_direct"
+ ):
raise NotImplementedError(
"eri_contraction='direct' is currently only available for PySCF."
)
@@ -254,16 +276,18 @@ def nuclear_gradient(excitation_or_mp, conv_tol=1e-9, eri_contraction=None,
with timer.record("contract_integral_derivatives"):
if eri_contraction == "direct":
grad = provider.correlated_gradient_direct(
- g1_ao, w_ao, g2_total, refstate=hf,
+ g1_ao,
+ w_ao,
+ g2_total,
+ refstate=hf,
shell_chunk_size=eri_shell_chunk_size,
pair_chunk_size=eri_pair_chunk_size,
pair_density_storage=eri_pair_density_storage,
)
else:
- grad = provider.correlated_gradient(
- g1_ao, w_ao, g2_ao_1, g2_ao_2
- )
+ grad = provider.correlated_gradient(g1_ao, w_ao, g2_ao_1, g2_ao_2)
- ret = GradientResult(excitation_or_mp, grad, g1, g2_total,
- timer, g1a=g1a, g2a=g2a)
+ ret = GradientResult(
+ excitation_or_mp, grad, g1, g2_total, timer, g1a=g1a, g2a=g2a
+ )
return ret
diff --git a/adcc/tests/geomopt_scanner_test.py b/adcc/tests/geomopt_scanner_test.py
index 0806b697..e78eac26 100644
--- a/adcc/tests/geomopt_scanner_test.py
+++ b/adcc/tests/geomopt_scanner_test.py
@@ -43,14 +43,18 @@
def _h2o_scf():
from pyscf import gto, scf
+
mol = gto.M(
atom="""
O 0 0 0
H 0 0 1.795239827225189
H 1.693194615993441 0 -0.599043184453037
""",
- basis="sto-3g", unit="Bohr", symmetry=False,
- verbose=0, parse_arg=False,
+ basis="sto-3g",
+ unit="Bohr",
+ symmetry=False,
+ verbose=0,
+ parse_arg=False,
)
mf = scf.RHF(mol)
mf.conv_tol = 1e-11
@@ -79,7 +83,9 @@ def test_scanner_validates_coordinate_shape():
def test_ground_state_scanner_matches_explicit_gradient_loop():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="mp2", gradient_kwargs={"eri_contraction": "full_ao"},
+ _h2o_scf(),
+ method="mp2",
+ gradient_kwargs={"eri_contraction": "full_ao"},
)
coords = scanner.initial_coords
@@ -94,7 +100,9 @@ def test_ground_state_scanner_matches_explicit_gradient_loop():
def test_calc_new_returns_geometric_engine_shape():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="mp2", gradient_kwargs={"eri_contraction": "full_ao"},
+ _h2o_scf(),
+ method="mp2",
+ gradient_kwargs={"eri_contraction": "full_ao"},
)
result = scanner.calc_new(scanner.initial_coords.ravel())
assert set(result) == {"energy", "gradient"}
@@ -107,7 +115,9 @@ def test_scanner_accepts_pyscf_mole_directly():
# passed straight to as_pyscf_method without a hand-rolled wrapper.
scfres = _h2o_scf()
scanner = adcc.NuclearGradientScanner(
- scfres, method="mp2", gradient_kwargs={"eri_contraction": "full_ao"},
+ scfres,
+ method="mp2",
+ gradient_kwargs={"eri_contraction": "full_ao"},
)
e_via_mole, g_via_mole = scanner(scfres.mol)
e_via_coords, g_via_coords = scanner(scanner.initial_coords)
@@ -117,10 +127,14 @@ def test_scanner_accepts_pyscf_mole_directly():
def test_scanner_step_callback_is_invoked_with_energy_and_gradient():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="mp2", gradient_kwargs={"eri_contraction": "full_ao"},
+ _h2o_scf(),
+ method="mp2",
+ gradient_kwargs={"eri_contraction": "full_ao"},
)
seen = []
- scanner.step_callback = lambda e, g: seen.append((e, float(np.asarray(g).sum())))
+ scanner.step_callback = lambda e, g: seen.append(
+ (e, float(np.asarray(g).sum()))
+ )
e, g = scanner(scanner.initial_coords)
assert len(seen) == 1
assert seen[0][0] == pytest.approx(e, abs=1e-12)
@@ -129,8 +143,12 @@ def test_scanner_step_callback_is_invoked_with_energy_and_gradient():
def test_run_adc_kwargs_are_forwarded_with_native_names():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="adc2", n_singlets=2, conv_tol=1e-7,
- follow="index", gradient_kwargs={"eri_contraction": "full_ao"},
+ _h2o_scf(),
+ method="adc2",
+ n_singlets=2,
+ conv_tol=1e-7,
+ follow="index",
+ gradient_kwargs={"eri_contraction": "full_ao"},
)
assert scanner.target.kwargs()["conv_tol"] == pytest.approx(1e-7)
assert scanner.target.kwargs()["n_singlets"] == 2
@@ -158,6 +176,7 @@ class FakeExcitation:
def test_unconverged_scf_raises():
from pyscf import scf
+
mf = scf.RHF(_h2o_scf().mol)
mf.max_cycle = 1
mf.conv_tol = 1e-12
@@ -169,8 +188,12 @@ def test_unconverged_scf_raises():
def test_excited_state_scanner_matches_explicit_gradient_loop():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="adc2", n_singlets=3, state_index=0,
- follow="index", conv_tol=1e-9,
+ _h2o_scf(),
+ method="adc2",
+ n_singlets=3,
+ state_index=0,
+ follow="index",
+ conv_tol=1e-9,
gradient_kwargs={"eri_contraction": "full_ao"},
)
@@ -180,10 +203,12 @@ def test_excited_state_scanner_matches_explicit_gradient_loop():
assert scanner.last_excitation.index == 0
assert scanner.last_tracking is None
- states = adcc.run_adc(scanner.last_scf, method="adc2", n_singlets=3,
- conv_tol=1e-9)
- explicit = adcc.nuclear_gradient(states.excitations[0],
- eri_contraction="full_ao")
+ states = adcc.run_adc(
+ scanner.last_scf, method="adc2", n_singlets=3, conv_tol=1e-9
+ )
+ explicit = adcc.nuclear_gradient(
+ states.excitations[0], eri_contraction="full_ao"
+ )
assert energy == pytest.approx(states.excitations[0].total_energy, abs=1e-7)
assert gradient.shape == (3, 3)
assert_allclose(gradient, explicit.total, atol=1e-7)
@@ -191,8 +216,12 @@ def test_excited_state_scanner_matches_explicit_gradient_loop():
def test_overlap_tracking_follows_same_state_character():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="adc2", n_singlets=3, state_index=1,
- follow="overlap", conv_tol=1e-9,
+ _h2o_scf(),
+ method="adc2",
+ n_singlets=3,
+ state_index=1,
+ follow="overlap",
+ conv_tol=1e-9,
gradient_kwargs={"eri_contraction": "full_ao"},
)
@@ -213,8 +242,12 @@ def test_overlap_tracking_follows_same_state_character():
def test_overlap_tracking_continuous_under_small_displacement():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="adc2", n_singlets=3, state_index=0,
- follow="overlap", conv_tol=1e-9,
+ _h2o_scf(),
+ method="adc2",
+ n_singlets=3,
+ state_index=0,
+ follow="overlap",
+ conv_tol=1e-9,
gradient_kwargs={"eri_contraction": "full_ao"},
)
@@ -233,8 +266,13 @@ def test_overlap_tracking_continuous_under_small_displacement():
def test_overlap_tracking_below_min_score_raises():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="adc2", n_singlets=3, state_index=0,
- follow="overlap", tracking_min_score=2.0, conv_tol=1e-9,
+ _h2o_scf(),
+ method="adc2",
+ n_singlets=3,
+ state_index=0,
+ follow="overlap",
+ tracking_min_score=2.0,
+ conv_tol=1e-9,
gradient_kwargs={"eri_contraction": "full_ao"},
)
scanner(scanner.initial_coords) # seed
@@ -244,8 +282,13 @@ def test_overlap_tracking_below_min_score_raises():
def test_overlap_tracking_ambiguous_gap_raises():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="adc2", n_singlets=3, state_index=0,
- follow="overlap", tracking_min_gap=10.0, conv_tol=1e-9,
+ _h2o_scf(),
+ method="adc2",
+ n_singlets=3,
+ state_index=0,
+ follow="overlap",
+ tracking_min_gap=10.0,
+ conv_tol=1e-9,
gradient_kwargs={"eri_contraction": "full_ao"},
)
scanner(scanner.initial_coords) # seed
@@ -255,7 +298,9 @@ def test_overlap_tracking_ambiguous_gap_raises():
def test_scf_guess_continuity_updates_previous_state():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="mp2", gradient_kwargs={"eri_contraction": "full_ao"},
+ _h2o_scf(),
+ method="mp2",
+ gradient_kwargs={"eri_contraction": "full_ao"},
)
energy_first, _ = scanner(scanner.initial_coords)
@@ -275,6 +320,7 @@ def test_scf_guess_continuity_updates_previous_state():
# Root-reorder tracking (FINDING 4)
# ---------------------------------------------------------------------------
+
class _FakeAoOperator:
def __init__(self, matrix):
self._matrix = matrix
@@ -295,7 +341,10 @@ def test_overlap_tracking_follows_reordered_root_by_overlap():
# index differs from the originally requested ``state_index``: the seeded
# descriptor matches candidate #1, while the requested state_index is 0.
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="adc2", n_singlets=3, state_index=0,
+ _h2o_scf(),
+ method="adc2",
+ n_singlets=3,
+ state_index=0,
follow="overlap",
)
nao = scanner.base_mol.nao
@@ -310,7 +359,8 @@ def projector(i):
seed_diffdm = projector(1)
scanner.previous_descriptor = _TrackingDescriptor(
mol=scanner.base_mol.copy(),
- transition_dm=seed_transition, state_diffdm=seed_diffdm,
+ transition_dm=seed_transition,
+ state_diffdm=seed_diffdm,
)
scanner.previous_index = 0
@@ -345,11 +395,15 @@ def projector(i):
# Validation / construction-time branches (FINDING 10)
# ---------------------------------------------------------------------------
-@pytest.mark.parametrize("kwargs", [
- {"target": "mp3"},
- {"method": "mp3"},
- {"target": GroundStateTarget(level=3)},
-])
+
+@pytest.mark.parametrize(
+ "kwargs",
+ [
+ {"target": "mp3"},
+ {"method": "mp3"},
+ {"target": GroundStateTarget(level=3)},
+ ],
+)
def test_ground_state_level_other_than_two_is_rejected(kwargs):
with pytest.raises(NotImplementedError, match="MP2 ground-state"):
adcc.NuclearGradientScanner(_h2o_scf(), **kwargs)
@@ -358,14 +412,21 @@ def test_ground_state_level_other_than_two_is_rejected(kwargs):
def test_invalid_follow_raises_eagerly_at_construction():
with pytest.raises(ValueError, match="overlap.*index"):
adcc.NuclearGradientScanner(
- _h2o_scf(), method="adc2", n_singlets=2, follow="bogus",
+ _h2o_scf(),
+ method="adc2",
+ n_singlets=2,
+ follow="bogus",
)
def test_state_index_out_of_range_raises():
scanner = adcc.NuclearGradientScanner(
- _h2o_scf(), method="adc2", n_singlets=2, state_index=5,
- follow="index", conv_tol=1e-9,
+ _h2o_scf(),
+ method="adc2",
+ n_singlets=2,
+ state_index=5,
+ follow="index",
+ conv_tol=1e-9,
gradient_kwargs={"eri_contraction": "full_ao"},
)
with pytest.raises(ValueError, match="out of range"):
@@ -379,10 +440,13 @@ def test_empty_excitations_raises_runtime_error():
scanner._select_excitation(states, scanner.base_mol.copy())
-@pytest.mark.parametrize("kwargs", [
- {"method": "mp2", "n_singlets": 3},
- {"target": "mp2", "n_singlets": 3},
-])
+@pytest.mark.parametrize(
+ "kwargs",
+ [
+ {"method": "mp2", "n_singlets": 3},
+ {"target": "mp2", "n_singlets": 3},
+ ],
+)
def test_ground_state_run_adc_kwargs_emit_runtime_warning(kwargs):
with pytest.warns(RuntimeWarning, match="Ignoring run_adc"):
adcc.NuclearGradientScanner(_h2o_scf(), **kwargs)
@@ -392,6 +456,7 @@ def test_ground_state_run_adc_kwargs_emit_runtime_warning(kwargs):
# Density-overlap scoring branches (FINDING 10)
# ---------------------------------------------------------------------------
+
def test_safe_ao_ndarray_swallows_unavailable_channels():
class FailingOperator:
def to_ndarray(self):
@@ -412,10 +477,12 @@ def state_diffdm_ao(self):
def test_tracking_score_records_unavailable_and_raises_when_both_missing():
scanner = adcc.NuclearGradientScanner(_h2o_scf(), method="adc2")
mol = scanner.base_mol.copy()
- previous = _TrackingDescriptor(mol=mol.copy(),
- transition_dm=None, state_diffdm=None)
- current = _TrackingDescriptor(mol=mol.copy(),
- transition_dm=None, state_diffdm=None)
+ previous = _TrackingDescriptor(
+ mol=mol.copy(), transition_dm=None, state_diffdm=None
+ )
+ current = _TrackingDescriptor(
+ mol=mol.copy(), transition_dm=None, state_diffdm=None
+ )
unavailable = set()
with pytest.raises(RuntimeError, match="neither transition_dm_ao"):
scanner._tracking_score(previous, current, mol, unavailable)