Relaxed coordinate scans

vq.relaxed_scan(...) drives one internal coordinate — a bond, an angle, or a dihedral — across a sequence of target values and relaxes all other degrees of freedom at each step. This is the canonical “scan the reaction coordinate” workflow used to locate transition-state guesses and map reaction energy profiles. It works for both molecules and periodic systems with the same API.

Quick start — molecular O–H stretch

import numpy as np
import vibeqc as vq

h2o = vq.Molecule(
    [
        vq.Atom(8, [0.00,  0.00, 0.0]),
        vq.Atom(1, [1.81,  0.00, 0.0]),
        vq.Atom(1, [-0.46, 1.75, 0.0]),
    ],
    charge=0, multiplicity=1,
)

result = vq.relaxed_scan(
    h2o, "sto-3g",
    coordinate=("bond", 0, 1),          # O (atom 0) → H (atom 1)
    values=np.linspace(1.6, 2.4, 9),    # bohr
    method="RHF",
    output="oh_stretch",                # writes oh_stretch.qvf
)

print(result.values)        # target bond lengths
print(result.energies)      # relaxed energies (Ha)
print(result.converged_flags)  # bool per scan point
print(result.geometries[0]) # vq.Molecule at the first scan point

The optimizer freezes the two constraint atoms (so the bond length stays exactly at the target) and lets every other atom relax. Each scan point warm-starts from the previous relaxed geometry, which gives smoother paths.

Quick start — periodic bond scan

The same call signature works for periodic systems; pass a PeriodicSystem and a kpoints mesh.

import numpy as np
import vibeqc as vq

L = np.eye(3) * 12.0
h2_box = vq.PeriodicSystem(
    3, L,
    [vq.Atom(1, [6.0, 6.0, 6.0]),
     vq.Atom(1, [7.4, 6.0, 6.0])],
    charge=0, multiplicity=1,
)

result = vq.relaxed_scan(
    h2_box, "sto-3g",
    coordinate=("bond", 0, 1),
    values=np.linspace(1.2, 2.4, 7),
    method="RHF",
    kpoints=(1, 1, 1),     # Γ-only
    output="h2_dissociation",
)

Coordinate types

Three coordinate kinds are supported. Bond targets are in bohr; angle and dihedral targets are in radians.

Coordinate	Form	Unit
Bond length	("bond", i, j)	bohr
Bond angle (vertex j)	("angle", i, j, k)	radians
Dihedral (i–j–k–l torsion)	("dihedral", i, j, k, l)	radians

Atom indices are 0-based and refer to the order of atoms in molecule.atoms (or system.unit_cell for periodic systems).

Freezing additional atoms — the slab pattern

For surface catalysis, freeze the slab bottom while the adsorbate moves. The constraint atoms are always frozen; pass extra indices via freeze_indices:

slab_info = vq.slab(...)     # see slabs_and_adsorbates.md
n_layers_frozen = 2

result = vq.relaxed_scan(
    slab_info.system, "def2-svp",
    coordinate=("bond", n2_idx_a, n2_idx_b),  # N2 bond stretch
    values=np.linspace(2.1, 4.0, 12),
    method="RKS", functional="PBE",
    kpoints=(2, 2, 1),
    freeze_indices=slab_info.bottom_layer_indices(n_layers_frozen),
    output="n2_dissociation_fe100",
)

The frozen-atom forwarding mirrors vibeqc.relax_atoms() (see CRYSTAL-style BIPOLE periodic HF and DFT).

Approach: rigid step with frozen endpoints

At each scan point:

1. The starting geometry is mutated so the constrained coordinate hits the target value exactly. For a bond, atom j is moved along the i → j direction. For an angle, atom k is rotated about vertex j. For a dihedral, atom l is rotated about the j → k axis.

2. The atoms participating in the constraint are added to the optimizer’s freeze_indices. The remaining atoms relax around the constrained geometry using the existing analytic-gradient L-BFGS-B driver.

3. The resulting geometry, energy, and convergence flag are recorded.

This is the same “rigid step” approach used by most QC packages for relaxed scans. The slightly more elegant alternative — projecting the gradient orthogonal to the constraint direction — is a planned v1.0 enhancement.

Output: ScanResult and the QVF reaction path

relaxed_scan returns a vibeqc.ScanResult:

result.values             # np.ndarray of target coordinate values
result.energies           # np.ndarray of relaxed energies (Ha)
result.geometries         # list of Molecule / PeriodicSystem
result.converged_flags    # list of bool
result.coordinate         # the constraint spec passed in
result.output_path        # Path to the QVF, if output= was set

When output="..." is set (or you call result.write_qvf(path)), a vibe-view-compatible reaction.path QVF is emitted. Frames become the scan animation; the reaction coordinate axis is the constraint value; endpoints are tagged as reactant / product waypoints, and the highest-energy frame is tagged as a generic point waypoint. (It is not tagged as a transition state — that requires a Hessian check; see hessians.)

Open the resulting .qvf with vibe-view to scrub the scan interactively.

Periodic scans → QVF v2. When geometries is a list of PeriodicSystem (periodic scan), the archive ships as qvf_version: 2 with the per-frame lattice + dim attached to the reaction.path section. vibe-view’s renderer (QVF Inc B) draws the unit-cell wireframe and wraps atom positions across in-plane periodic boundaries — slab scans animate with the cell visible and adsorbates don’t jump across the cell boundary mid-scan.

Citation surface. write_qvf emits a citations section containing the BibTeX for every per-frame SCF (method / basis / functional). The relaxed-scan driver itself doesn’t fire a driver-specific paper (relaxed coordinate scans are textbook; there is no single “scan paper” to cite). NEB callers needing the Henkelman+Jónsson 2000 / Henkelman+Uberuaga+Jónsson 2000 / Smidstrup 2014 citations should use vibeqc.NEBResult.write_qvf() instead — it fires the routes.drivers.neb (and optionally ci_neb) routes automatically.

Limits and out-of-scope features

• 2D scans (sweeping two coordinates simultaneously) are not yet supported. Run two nested 1D scans by hand for now.

• Constrained-gradient projection is not implemented; the rigid- step approach above is what ships today.

• Synchronous-transit / NEB path methods live in a separate workstream — relaxed_scan is the simpler reaction-coordinate driver.