Construct a Reaction-Path Problem

A reaction-path problem bundles all computational ingredients needed to evaluate the quantum-mechanical ground-state energy of a sequence of molecular geometries (reaction path or scan). Each geometry is evaluated as a separate ground-state problem, with optional specialized handling for embedded calculations.

Goal

Build a reaction-path problem from a reaction configuration.

The reaction-path problem is a key component of a typical workflow for a quantum chemistry calculation which follows three main steps:

1. Define the molecular configuration: Specify the atoms, their coordinates, and the basis set. This creates a molecular configuration that encodes the physical system you want to study.

2. Build the problem: Use the molecular configuration to construct a problem object (such as a ground state energy problem). This bundles all necessary quantum-chemical data (integrals, electron counts, etc.) for the calculation.

3. Choose and run the calculator: Select a calculator (e.g., FAST-VQE, BEAST-VQE, ADAPT-VQE) that knows how to solve the problem. The calculator takes the problem object and computes the desired properties, such as the total energy.

This “molecule → problem → calculator” flow ensures a clear separation between system definition, problem formulation, and computational method, making it easy to swap components.

Prerequisites

• A reaction configuration (see Create a ReactionConfiguration)

• (For simple reaction builder) a ground-state problem builder (see Construct a Ground State Energy Problem)

Steps

1. Create the reaction-path problem builder (simple case)

   The reaction_path() builder wraps a ground-state problem builder to run it over all geometries in the reaction configuration.

   • .simple(problem_builder): applies the same ground-state problem builder to every geometry.

   See Understanding and Using Kvantify Qrunch’s Fluent Builder Pattern for details on the builder pattern and how to use them.

   Example:

   import qrunch as qc
   
   # Create a simple reaction-path builder
   reaction_problem_builder = (
       qc.problem_builder_creator()
       .reaction_path()
       .simple(problem_builder)  # Required ground-state problem builder
       .create()
   )
   

2. Build the reaction-path problem

   You can build either a restricted or unrestricted reaction-path problem, just like for single-geometry ground-state problems:

   # Unrestricted
   reaction_problem = reaction_problem_builder.build_unrestricted(reaction_configuration)
   
   # or Restricted
   reaction_problem = reaction_problem_builder.build_restricted(reaction_configuration)
   

   In both cases, the result is a reaction-path problem object containing one ground-state problem per geometry.

Alternative: Even-Handed Reaction Builder

The even_handed() reaction-path builder implements the Even-handed subsystem selection method [WMM18].

This approach:

• Uses Wavefunction-in-DFT projection-based embedding

• Selects subsystems consistently along the path

• Accepts all configuration options of the projective_embedding ground-state problem builder (see Create a Projective-Embedding Ground-State Problem for details)

Example:

import qrunch as qc

even_handed_reaction_problem_builder = (
    qc.problem_builder_creator()
    .reaction_path()
    .even_handed()
    # Add .<projective_embedding options> here
    .create()
)

reaction_problem = even_handed_reaction_problem_builder.build_restricted(reaction_configuration)

Consistent Active Space Along the Path

Selecting an active space independently at every geometry can pick different orbitals at different points along the path. The resulting energy curve then contains discontinuities that are artefacts of the active-space selection rather than of the chemistry. with_consistent_active_space() avoids this by determining one active-space definition from all configurations at once, and then applying it consistently to every configuration.

It is available on both reaction-path builders:

import qrunch as qc

problem_builder = (
    qc.problem_builder_creator()
    .ground_state()
    .standard()
    .choose_molecular_orbital_calculator()
    .moller_plesset_2()  # or .ccsd()
    .create()
)

reaction_problem_builder = (
    qc.problem_builder_creator()
    .reaction_path()
    .simple(problem_builder)
    .with_consistent_active_space()
    .create()
)

reaction_problem = reaction_problem_builder.build_restricted(reaction_configuration)

For the even-handed builder the call is placed in the same way:

even_handed_reaction_problem_builder = (
    qc.problem_builder_creator()
    .reaction_path()
    .even_handed()
    .with_consistent_active_space()
    .create()
)

How It Works

The procedure runs in two passes over the configurations:

1. Discover which atomic shells matter. At each configuration, the orbitals whose natural occupation deviates from the closed-shell ideal values (see Define an Active Space (Complete Active Space)) are taken as candidates. Each candidate is fingerprinted using its overlap with localized intrinsic bonding orbitals (IBOs) [Kni13], which labels it with the atomic shells it is built from, for example C 2p or N 2p. The union of the fingerprints from all configurations forms the master label set.

2. Build the same active space everywhere. At each configuration, AVAS (atomic valence active space) [SSCK17] projects the molecular orbitals onto the master label set, and the orbitals whose projection exceeds avas_threshold become active. Because every configuration is projected onto the same label set, the active space describes the same chemistry all along the path, even though the individual orbitals change shape as the geometry changes.

Requirements

• The ground-state problem builder must use a correlated molecular orbital calculator (MP2 or CCSD). A mean-field calculator only produces integer occupations, so no candidate orbitals can be identified, and an error is raised.

• For unrestricted problems, and when use_unrestricted_natural_orbitals=True, the molecular orbital calculator must be configured with spin_summed_natural_orbitals=True. See Common Options for the available options and a broken-symmetry example.

• You do not need to specify the size of the active space: it is determined by the procedure. If the underlying ground-state problem builder has an active_space(...) modifier, it is replaced by the consistent active space. Modifiers added before it are still applied before the active-space reduction, and modifiers added after it, such as to_dense_integrals(), are still applied afterwards.

Tuning the Consistent Active Space

You can pass a ConsistentActiveSpaceOptions instance to control the selection:

import qrunch as qc

options = qc.options.ConsistentActiveSpaceOptions(
    occupation_deviation_threshold=0.02,
    max_active_spatial_orbitals=12,
)

reaction_problem_builder = (
    qc.problem_builder_creator()
    .reaction_path()
    .simple(problem_builder)
    .with_consistent_active_space(options)
    .create()
)

Common Options

These are the options you will normally tune. They control which orbitals are considered and how large the resulting active space is allowed to become.

Option	Default	Description
occupation_deviation_threshold	0.02	Minimum occupation deviation for an orbital to become a candidate in pass 1. Raising it keeps only the most strongly partially occupied orbitals and gives a smaller active space.
max_active_spatial_orbitals	None	Hard cap on the number of active spatial orbitals per configuration. When the cap is exceeded, a stricter AVAS threshold is found by bisection and re-applied to all configurations, so the active space stays consistent. Mutually exclusive with num_active_spatial_orbitals.
num_active_spatial_orbitals	None	Ask for an active space of exactly this many spatial orbitals at every configuration. Mutually exclusive with max_active_spatial_orbitals: use that one if an upper bound is enough, and this one if you need a specific size, for example to compare several reaction paths on an equal footing or to fit a fixed qubit budget. The builder tunes the other options for you until it reaches the requested size, so you do not need to adjust them by hand. An exact hit is not always possible; if none is found, a warning tells you which size was used instead.
reference_basis	"ano-rcc-mb"	Atomic reference basis defining the IAOs used for fingerprinting and the target space AVAS projects onto. The default is a minimal basis, so the reference space is the valence only. The next rungs, "ano-rcc-vdzp" and "ano-rcc-vtzp", add correlating shells such as \(3s\), \(3p\) and \(3d\) to the reference, which the fingerprints can then select and AVAS can then project onto.
conjugated_orbitals_only	False	Keep only candidates whose fingerprints consist exclusively of \(np\) shells on non-hydrogen atoms. Use this when the chemically relevant active space is the \(\pi\) subsystem.
use_unrestricted_natural_orbitals	False	Obtain the natural orbitals for a restricted problem from an unrestricted calculation, so that options that only take effect in UHF/UMP2 can influence the occupation numbers. Setting this alone does not improve the description of a strongly correlated system: an unrestricted calculation started from a closed-shell guess simply converges back to the restricted solution. To get the broken-symmetry solution you normally also need break_spin_symmetry=True and/or stability_analysis=True on the molecular orbital calculator, together with spin_summed_natural_orbitals=True.

A broken-symmetry setup for a restricted problem therefore looks like this:

import qrunch as qc

problem_builder = (
    qc.problem_builder_creator()
    .ground_state()
    .standard()
    .choose_molecular_orbital_calculator()
    .moller_plesset_2(
        qc.options.MollerPlesset2CalculatorOptions(
            spin_summed_natural_orbitals=True,
            break_spin_symmetry=True,
            stability_analysis=True,
        )
    )
    .create()
)

reaction_problem_builder = (
    qc.problem_builder_creator()
    .reaction_path()
    .simple(problem_builder)
    .with_consistent_active_space(
        qc.options.ConsistentActiveSpaceOptions(use_unrestricted_natural_orbitals=True)
    )
    .create()
)

restricted_problem = reaction_problem_builder.build_restricted(reaction_configuration)

See MollerPlesset2CalculatorOptions and CCSDCalculatorOptions for the full list of molecular orbital calculator options.

Asking for a Specific Active-Space Size

By default the procedure picks the active-space size for you. If you need a specific size instead, set num_active_spatial_orbitals:

import qrunch as qc

reaction_problem_builder = (
    qc.problem_builder_creator()
    .reaction_path()
    .simple(problem_builder)
    .with_consistent_active_space(
        qc.options.ConsistentActiveSpaceOptions(num_active_spatial_orbitals=10)
    )
    .create()
)

The builder then tries different settings on your behalf until every configuration along the path comes out with ten active spatial orbitals. It starts from the values you configured and adjusts them automatically, so in most cases this is the only option you need to set. avas_threshold is ignored while this is active, because the builder determines it itself.

reference_basis is one of the settings it varies. If no combination of thresholds reaches the requested size against the configured reference, the builder climbs the nested ANO-RCC ladder, ano-rcc-mb then ano-rcc-vdzp then ano-rcc-vtzp. Each rung adds correlating shells that the fingerprints can select, which raises the largest active space that can be reached at all. A rung that does not suit the molecule is skipped and reported at WARNING level; ano-rcc-vtzp in particular is only usable when the orbital basis is at least triple zeta, since otherwise it contributes more reference orbitals than the calculation has.

The size can only change in whole orbitals, and sometimes it jumps past the number you asked for, so an exact hit is not guaranteed. When that happens the builder keeps the closest active space that is still the same at every configuration and logs a warning telling you which size it used and which settings produced it, including the reference basis. If the size it reports is not acceptable, ask for a nearby number instead.

Expert Options

These options control the internals of the fingerprinting and projection steps, and rarely need changing.

Option	Default	Description
population_contribution_threshold	0.1	Minimum fraction of a candidate’s density that must sit on an atom-shell for that label to be kept. Raising it gives sparser fingerprints and a smaller master label set.
overlap_contribution_threshold	0.1	Minimum squared overlap between a candidate orbital and an IBO for that IBO to contribute to the fingerprint. Lower it for delocalized \(\pi\) systems where no single IBO dominates.
avas_threshold	0.2	Minimum projection onto the master label set for an orbital to be made active in pass 2. Lowering it enlarges the active space and reduces the risk of different configurations picking different numbers of active orbitals.

Verify the Result

• At INFO log level, Qrunch reports the master label set and, for each configuration, how many active orbitals and electrons AVAS selected. These numbers should be the same for all configurations.

• At DEBUG log level, Qrunch also reports the per-configuration fingerprint tables, which show why a given atom-shell label ended up in the master label set.

See Also

• Create a ReactionConfiguration

• Construct a Ground State Energy Problem

• Define an Active Space (Complete Active Space)

• Create a Projective-Embedding Ground-State Problem

• API reference: reaction_path_problem

Next Step

You can use the reaction-path problem to calculate energies along the path:

• Calculate the Ground State Energy Using the FAST-VQE

---

References

[Kni13, SSCK17, WMM18]