Performing a Reaction-Path Potential Energy Surface (PES) study
In this example, we build a Potential Energy Surface (PES) for a localized dissociation — the O–H stretch in methanol — while keeping the methyl group as a spectator. We compare three ways to assemble a reaction path in Qrunch:
Simple + Standard — apply the same ground-state builder independently to each geometry.
Simple + Projective-Embedding — still independent by-geometry, but using a Projective-Embedding ground-state builder.
Even-Handed Projective-Embedding — a reaction-path-aware Projective-Embedding strategy that enforces consistent embedded subsystems across all geometries (see https://arxiv.org/abs/1809.03004) to avoid discontinuities/cusps in the PES.
The broader motivation for projection-based embedding is described in Projection-Based Wavefunction-in-DFT Embedding.
Full script
The script below shows the full script, which are then explained section by section.
1"""Demonstrate calculation of a Potential Energy Surface (PES) of a dissociation reaction."""
2# pyright: reportUnknownVariableType=false, reportMissingImports=false, reportCallIssue=false, reportUnknownArgumentType=false
3# ruff: noqa
4
5from typing import Literal, TypeAlias
6
7from pathlib import Path
8
9import matplotlib.pyplot as plt
10import numpy as np
11from numpy.typing import NDArray
12import qrunch as qc
13
14
15qc.setup_logger(qc.LOGGER_INFO)
16AtomSymbol = Literal["C", "O", "H"]
17Atom: TypeAlias = tuple[AtomSymbol, float, float, float]
18
19
20def methanol_with_stretched_bond(new_distance: float) -> list[Atom]:
21 """
22 Provide a methanol molecule with a stretched OH bond, by translating the Hydrogen along the bond direction.
23
24 Args:
25 new_distance: Target O-H distance in Å.
26
27 """
28 atoms: list[Atom] = [
29 ("C", 0.0000, 0.0000, 0.0000),
30 ("O", 1.4300, 0.0000, 0.0000), # index 1 (O)
31 ("H", 1.9300, 0.9300, 0.0000), # index 2 (OH hydrogen) <-- this one moves
32 ("H", -0.5400, 0.9400, 0.0000),
33 ("H", -0.5400, -0.4700, 0.8150),
34 ("H", -0.5400, -0.4700, -0.8150),
35 ]
36 oxygen_index = 1
37 hydrogen_index = 2
38
39 _, x_h, y_h, z_h = atoms[hydrogen_index]
40 _, x_o, y_o, z_o = atoms[oxygen_index]
41
42 vector = np.array([x_h - x_o, y_h - y_o, z_h - z_o], dtype=float)
43 vector_norm = float(np.linalg.norm(vector))
44 direction_vector = vector / vector_norm
45 new_position = np.array([x_o, y_o, z_o], dtype=float) + new_distance * direction_vector
46
47 stretched_atoms = list(atoms)
48 new_h: Atom = ("H", float(new_position[0]), float(new_position[1]), float(new_position[2]))
49 stretched_atoms[hydrogen_index] = new_h
50 return stretched_atoms
51
52
53
54
55def plot_relative_energy_vs_coordinate(
56 reaction_coordinates: NDArray[np.float64],
57 energies_hartree: list[float],
58 pe_energies_hartree: list[float],
59 even_handed_energies_hartree: list[float],
60 *,
61 outdir: Path = Path("dist"),
62 show: bool = False,
63) -> None:
64 """
65 Plot relative energy (kcal/mol) vs reaction coordinate and save PNG and CSV.
66
67 The relative energy is referenced to the minimum energy among the provided points.
68
69 Args:
70 reaction_coordinates: Reaction coordinates (Å).
71 energies_hartree: Total energies (Hartree) corresponding to q_values.
72 pe_energies_hartree: Total energies (Hartree) corresponding to q_values for Projected embedding.
73 even_handed_energies_hartree: Total energies (Hartree) corresponding to q_values for even handed.
74 outdir: Output directory where plots are written. Defaults to "dist".
75 show: Whether to display the figures interactively after saving.
76
77 """
78 outdir.mkdir(parents=True, exist_ok=True)
79
80 oh_distances = np.asarray(list(reaction_coordinates), dtype=float)
81 energies = np.asarray(list(energies_hartree), dtype=float)
82 pe_energies = np.asarray(list(pe_energies_hartree), dtype=float)
83 even_handed_energies = np.asarray(list(even_handed_energies_hartree), dtype=float)
84
85 # Compute relative energies (subtract min)
86 rel_simple = energies - np.min(energies)
87 rel_pe = pe_energies - np.min(pe_energies)
88 rel_even = even_handed_energies - np.min(even_handed_energies)
89
90 plt.figure(figsize=(6.0, 4.0))
91 plt.plot(oh_distances, rel_simple, "-o", label="Simple+Standard", color="darkgreen")
92 plt.plot(oh_distances, rel_pe, "--o", label="Simple+Projected-Embedding", color="green")
93 plt.plot(oh_distances, rel_even, ":o", label="Even-Handed", color="lightgreen")
94 plt.xlabel(r"O-H distance [Å]")
95 plt.ylabel("Energy [Hartree]")
96 plt.title("Methanol dissociation reaction PES")
97 plt.grid(visible=True, alpha=0.3)
98 plt.legend()
99 plt.tight_layout()
100
101 png_path = outdir / "methanol_rxn.png"
102 plt.savefig(png_path, dpi=200)
103 if show:
104 plt.show()
105 else:
106 plt.close()
107
108
109def main(oh_distances: NDArray[np.float64], data_path: Path | None = None) -> list[float]:
110 """
111 Scan methanol O-H dissociation where only the hydroxyl H moves.
112
113 Args:
114 oh_distances: Array of O-H distances [Å] to evaluate.
115
116 Returns:
117 List of total energies [Hartree] corresponding to the input distances.
118
119 """
120 if data_path is None:
121 directory = Path.cwd() / "data"
122 else:
123 directory = data_path
124
125 # =======================================================================
126 # Simple Reaction Path Problem Builder with standard ground state problem
127 # =======================================================================
128
129 ground_state_problem_builder = (
130 qc.problem_builder_creator() # Start creating a problem builder
131 .ground_state() # Narrow to a ground state problem
132 .standard() # Narrow to a standard ground state problem
133 .choose_molecular_orbital_calculator()
134 .hartree_fock(
135 options=qc.options.HartreeFockCalculatorOptions(do_density_fitting=True)
136 ) # Use Hartree-Fock to construct molecular orbitals
137 .choose_repulsion_integral_builder()
138 .resolution_of_the_identity(auxiliary_basis="weigend")
139 .add_problem_modifier()
140 .active_space( # Restrict to an active space (16,8)
141 number_of_active_spatial_orbitals=8,
142 number_of_active_alpha_electrons=4,
143 )
144 .add_problem_modifier()
145 .to_dense_integrals()
146 .choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
147 .file_persister( # Perform the sub-choice selection - Here we pick file persister
148 directory=directory, # Directory where data files will be saved
149 extension=".qdk", # File extension for the data files
150 do_save=True, # Whether to save data to files
151 load_policy="fallback", # Load data if available, otherwise compute
152 overwriting_policy="rename", # Rename files if filename already exist
153 padding_width=3, # Padding width for file numbering (001, 002, ...)
154 )
155 .create()
156 )
157
158 # Build the reaction path problem builder.
159 reaction_problem_builder = (
160 qc.problem_builder_creator() # Start creating a problem builder
161 .reaction_path() # Narrow: pick ground state problem
162 .simple( # Narrow: pick simple reaction path problem
163 ground_state_problem_builder # Use the ground state problem builder for each geometry
164 )
165 .create()
166 )
167
168 # ===================================================================================
169 # Simple Reaction Path Problem Builder with projective embedding ground state problem
170 # ===================================================================================
171
172 # Build Projective embedding ground state problem builder.
173 pe_ground_state_problem_builder = (
174 qc.problem_builder_creator() # Start creating a problem builder
175 .ground_state() # Narrow to a ground state problem
176 .projective_embedding() # Narrow to a Projective embedding ground state problem
177 .choose_embedded_orbital_calculator()
178 .hartree_fock(
179 options=qc.options.HartreeFockCalculatorOptions(do_density_fitting=True)
180 ) # Use Hartree-Fock to construct molecular orbitals
181 .choose_repulsion_integral_builder()
182 .resolution_of_the_identity(auxiliary_basis="weigend")
183 .choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
184 .file_persister( # Perform the sub-choice selection - Here we pick file persister
185 directory=directory, # Directory where data files will be saved
186 extension=".qdk", # File extension for the data files
187 do_save=True, # Whether to save data to files
188 load_policy="fallback", # Load data if available, otherwise compute
189 overwriting_policy="rename", # Rename files if filename already exist
190 padding_width=3, # Padding width for file numbering (001, 002, ...)
191 )
192 .add_problem_modifier()
193 .to_dense_integrals()
194 .create()
195 )
196
197 # Build the reaction path problem builder
198 # with the Projective embedding ground state problem builder.
199 pe_reaction_problem_builder = (
200 qc.problem_builder_creator() # Start creating a problem builder
201 .reaction_path() # Narrow: pick ground state problem
202 .simple( # Narrow: pick simple reaction path problem
203 pe_ground_state_problem_builder
204 )
205 .create()
206 )
207
208 # ===========================================================================
209 # Even handed Reaction Path Problem Builder https://arxiv.org/abs/1809.03004
210 # ===========================================================================
211
212 # Build the reaction problem builder.
213 even_handed_reaction_problem_builder = (
214 qc.problem_builder_creator() # Start creating a problem builder
215 .reaction_path() # Narrow: pick ground state problem
216 .even_handed() # Narrow: pick even-handed reaction path problem
217 .choose_embedded_orbital_calculator()
218 .hartree_fock(
219 options=qc.options.HartreeFockCalculatorOptions(do_density_fitting=True)
220 ) # Use Hartree-Fock to construct molecular orbitals
221 .choose_repulsion_integral_builder()
222 .resolution_of_the_identity(auxiliary_basis="weigend")
223 .choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
224 .file_persister( # Perform the sub-choice selection - Here we pick file persister
225 directory=directory, # Directory where data files will be saved
226 extension=".qdk", # File extension for the data files
227 do_save=True, # Whether to save data to files
228 load_policy="fallback", # Load data if available, otherwise compute
229 overwriting_policy="rename", # Rename files if filename already exist
230 padding_width=3, # Padding width for file numbering (001, 002, ...)
231 )
232 .add_problem_modifier()
233 .to_dense_integrals()
234 .create()
235 )
236
237 # Build the FAST-VQE calculator
238 fast_vqe_calculator = (
239 qc.calculator_creator() # Start creating a calculator
240 .vqe() # Narrow: pick Variational quantum eigensolver (VQE)
241 .iterative() # Narrow: pick the iterative VQE
242 .standard() # Narrow: pick the standard ansatz VQE (FAST-VQE)
243 .with_options(options=qc.options.IterativeVqeOptions(max_iterations=200))
244 .choose_stopping_criterion()
245 .patience(patience=2, threshold=1e-3)
246 .choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
247 .file_persister( # Perform the sub-choice selection - Here we pick file persister
248 directory=directory, # Directory where data files will be saved
249 extension=".qdk", # File extension for the data files
250 do_save=True, # Whether to save data to files
251 load_policy="fallback", # Load data if available, otherwise compute
252 overwriting_policy="rename", # Rename files if filename already exist
253 padding_width=3, # Padding width for file numbering (001, 002, ...)
254 )
255 .create() # Create the calculator instance
256 )
257
258 # Build the list of molecules that make up the reaction path.
259 reaction = []
260 for distance in oh_distances:
261 molecule = methanol_with_stretched_bond(distance)
262 reaction.append(molecule)
263
264 reaction_configuration = qc.build_reaction_configuration(
265 reaction=reaction,
266 basis_set="sto3g",
267 spin_difference=0,
268 charge=0,
269 embedded_atoms=[1, 2],
270 aux_basis_set="weigend",
271 )
272
273 reaction_path_problem = reaction_problem_builder.build_unrestricted(reaction_configuration)
274 reaction_result = fast_vqe_calculator.calculate(reaction_path_problem)
275
276 pe_reaction_path_problem = pe_reaction_problem_builder.build_unrestricted(reaction_configuration)
277 pe_reaction_result = fast_vqe_calculator.calculate(pe_reaction_path_problem)
278
279 even_handed_reaction_problem = even_handed_reaction_problem_builder.build_unrestricted(reaction_configuration)
280 even_handed_reaction_result = fast_vqe_calculator.calculate(even_handed_reaction_problem)
281
282 energies = reaction_result.total_energies.values
283 pe_energies = pe_reaction_result.total_energies.values
284 even_handed_energies = even_handed_reaction_result.total_energies.values
285
286 # Plot and save to dist/
287 outdir = Path("dist")
288 plot_relative_energy_vs_coordinate(
289 reaction_coordinates=oh_distances,
290 energies_hartree=energies,
291 pe_energies_hartree=pe_energies,
292 even_handed_energies_hartree=even_handed_energies,
293 outdir=outdir,
294 )
295
296
297 return reaction_result.total_energies.values
298
299
300if __name__ == "__main__":
301 # OH distances in Å
302 distances = np.linspace(1.0, 3.4, 6)
303 main(distances)
Explanation
Imports and setup
from typing import Literal, TypeAlias
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from numpy.typing import NDArray
import qrunch as qc
import qrunch as qc loads the public, stable API.
For most users this is the only import needed, and the one guaranteed to remain stable across versions.
Defining the localized reaction coordinate
def methanol_with_stretched_bond(new_distance: float) -> list[Atom]:
"""
Provide a methanol molecule with a stretched OH bond, by translating the Hydrogen along the bond direction.
Args:
new_distance: Target O-H distance in Å.
"""
atoms: list[Atom] = [
("C", 0.0000, 0.0000, 0.0000),
("O", 1.4300, 0.0000, 0.0000), # index 1 (O)
("H", 1.9300, 0.9300, 0.0000), # index 2 (OH hydrogen) <-- this one moves
("H", -0.5400, 0.9400, 0.0000),
("H", -0.5400, -0.4700, 0.8150),
("H", -0.5400, -0.4700, -0.8150),
]
oxygen_index = 1
hydrogen_index = 2
_, x_h, y_h, z_h = atoms[hydrogen_index]
_, x_o, y_o, z_o = atoms[oxygen_index]
vector = np.array([x_h - x_o, y_h - y_o, z_h - z_o], dtype=float)
vector_norm = float(np.linalg.norm(vector))
direction_vector = vector / vector_norm
new_position = np.array([x_o, y_o, z_o], dtype=float) + new_distance * direction_vector
stretched_atoms = list(atoms)
new_h: Atom = ("H", float(new_position[0]), float(new_position[1]), float(new_position[2]))
stretched_atoms[hydrogen_index] = new_h
return stretched_atoms
We create each geometry by translating only the hydroxyl hydrogen along the O–H bond direction to the requested distance. The methyl group and the C–O bond remain fixed — ensuring that only the O–H pair participates in the reaction.
Simple reaction path: standard ground-state builder per geometry
# =======================================================================
# Simple Reaction Path Problem Builder with standard ground state problem
# =======================================================================
ground_state_problem_builder = (
qc.problem_builder_creator() # Start creating a problem builder
.ground_state() # Narrow to a ground state problem
.standard() # Narrow to a standard ground state problem
.choose_molecular_orbital_calculator()
.hartree_fock(
options=qc.options.HartreeFockCalculatorOptions(do_density_fitting=True)
) # Use Hartree-Fock to construct molecular orbitals
.choose_repulsion_integral_builder()
.resolution_of_the_identity(auxiliary_basis="weigend")
.add_problem_modifier()
.active_space( # Restrict to an active space (16,8)
number_of_active_spatial_orbitals=8,
number_of_active_alpha_electrons=4,
)
.add_problem_modifier()
.to_dense_integrals()
.choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
.file_persister( # Perform the sub-choice selection - Here we pick file persister
directory=directory, # Directory where data files will be saved
extension=".qdk", # File extension for the data files
do_save=True, # Whether to save data to files
load_policy="fallback", # Load data if available, otherwise compute
overwriting_policy="rename", # Rename files if filename already exist
padding_width=3, # Padding width for file numbering (001, 002, ...)
)
.create()
)
# Build the reaction path problem builder.
reaction_problem_builder = (
qc.problem_builder_creator() # Start creating a problem builder
.reaction_path() # Narrow: pick ground state problem
.simple( # Narrow: pick simple reaction path problem
ground_state_problem_builder # Use the ground state problem builder for each geometry
)
.create()
)
This Simple + Standard path applies the same ground-state problem to each geometry independently. To keep cost modest we:
request Hartree–Fock orbitals with density fitting (
do_density_fitting=True),use resolution-of-the-identity for repulsion integrals (
weigend),restrict to an active space (here 8 spatial orbitals with 4 alpha electrons),
store/reuse intermediates via a file persister.
Simple reaction path with a Projective-Embedding ground-state builder
# ===================================================================================
# Simple Reaction Path Problem Builder with projective embedding ground state problem
# ===================================================================================
# Build Projective embedding ground state problem builder.
pe_ground_state_problem_builder = (
qc.problem_builder_creator() # Start creating a problem builder
.ground_state() # Narrow to a ground state problem
.projective_embedding() # Narrow to a Projective embedding ground state problem
.choose_embedded_orbital_calculator()
.hartree_fock(
options=qc.options.HartreeFockCalculatorOptions(do_density_fitting=True)
) # Use Hartree-Fock to construct molecular orbitals
.choose_repulsion_integral_builder()
.resolution_of_the_identity(auxiliary_basis="weigend")
.choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
.file_persister( # Perform the sub-choice selection - Here we pick file persister
directory=directory, # Directory where data files will be saved
extension=".qdk", # File extension for the data files
do_save=True, # Whether to save data to files
load_policy="fallback", # Load data if available, otherwise compute
overwriting_policy="rename", # Rename files if filename already exist
padding_width=3, # Padding width for file numbering (001, 002, ...)
)
.add_problem_modifier()
.to_dense_integrals()
.create()
)
# Build the reaction path problem builder
# with the Projective embedding ground state problem builder.
pe_reaction_problem_builder = (
qc.problem_builder_creator() # Start creating a problem builder
.reaction_path() # Narrow: pick ground state problem
.simple( # Narrow: pick simple reaction path problem
pe_ground_state_problem_builder
)
.create()
)
This is Simple + Projective-Embedding: still geometry-by-geometry, but now the ground-state builder is the PE variant. In this minimal example we pick canonical Hartree–Fock orbitals from which to construct the Hamiltonian used by FAST-VQE. (All the usual PE options are available.)
Even-Handed Projective-Embedding reaction path
# ===========================================================================
# Even handed Reaction Path Problem Builder https://arxiv.org/abs/1809.03004
# ===========================================================================
# Build the reaction problem builder.
even_handed_reaction_problem_builder = (
qc.problem_builder_creator() # Start creating a problem builder
.reaction_path() # Narrow: pick ground state problem
.even_handed() # Narrow: pick even-handed reaction path problem
.choose_embedded_orbital_calculator()
.hartree_fock(
options=qc.options.HartreeFockCalculatorOptions(do_density_fitting=True)
) # Use Hartree-Fock to construct molecular orbitals
.choose_repulsion_integral_builder()
.resolution_of_the_identity(auxiliary_basis="weigend")
.choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
.file_persister( # Perform the sub-choice selection - Here we pick file persister
directory=directory, # Directory where data files will be saved
extension=".qdk", # File extension for the data files
do_save=True, # Whether to save data to files
load_policy="fallback", # Load data if available, otherwise compute
overwriting_policy="rename", # Rename files if filename already exist
padding_width=3, # Padding width for file numbering (001, 002, ...)
)
.add_problem_modifier()
.to_dense_integrals()
.create()
)
Even-Handed PE (https://arxiv.org/abs/1809.03004) constructs a consistent embedded subsystem across the entire reaction path by tracking localized orbitals and their partitioning jointly for the set of geometries. This prevents cusps/discontinuities that can occur when the embedded set changes abruptly along the path. It typically yields smoother, more physical PES curves.
Note
Every option available to the PE ground-state builder (e.g., embedded orbital calculator choices, localization and assignment choices, etc.) is also available in the Even-Handed reaction path builder.
FAST-VQE calculator
# Build the FAST-VQE calculator
fast_vqe_calculator = (
qc.calculator_creator() # Start creating a calculator
.vqe() # Narrow: pick Variational quantum eigensolver (VQE)
.iterative() # Narrow: pick the iterative VQE
.standard() # Narrow: pick the standard ansatz VQE (FAST-VQE)
.with_options(options=qc.options.IterativeVqeOptions(max_iterations=200))
.choose_stopping_criterion()
.patience(patience=2, threshold=1e-3)
.choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
.file_persister( # Perform the sub-choice selection - Here we pick file persister
directory=directory, # Directory where data files will be saved
extension=".qdk", # File extension for the data files
do_save=True, # Whether to save data to files
load_policy="fallback", # Load data if available, otherwise compute
overwriting_policy="rename", # Rename files if filename already exist
padding_width=3, # Padding width for file numbering (001, 002, ...)
)
.create() # Create the calculator instance
)
We use FAST-VQE (iterative) with a small patience and a loose threshold so the example runs quickly. For production, tighten these criteria. We also include a file persister to store intermediates.
Assemble the reaction and build a ReactionConfiguration
# Build the list of molecules that make up the reaction path.
reaction = []
for distance in oh_distances:
molecule = methanol_with_stretched_bond(distance)
reaction.append(molecule)
reaction_configuration = qc.build_reaction_configuration(
reaction=reaction,
basis_set="sto3g",
spin_difference=0,
charge=0,
embedded_atoms=[1, 2],
aux_basis_set="weigend",
)
We pass the list of geometries to qc.build_reaction_configuration(...),
set the basis, charge/spin, and — for PE/Even-Handed — declare the embedded atoms.
Here we embed the O and H of the hydroxyl group (indices 1 and 2) to focus correlation only where the bond is breaking.
Compute energies for each path
reaction_path_problem = reaction_problem_builder.build_unrestricted(reaction_configuration)
reaction_result = fast_vqe_calculator.calculate(reaction_path_problem)
pe_reaction_path_problem = pe_reaction_problem_builder.build_unrestricted(reaction_configuration)
pe_reaction_result = fast_vqe_calculator.calculate(pe_reaction_path_problem)
even_handed_reaction_problem = even_handed_reaction_problem_builder.build_unrestricted(reaction_configuration)
even_handed_reaction_result = fast_vqe_calculator.calculate(even_handed_reaction_problem)
This evaluates the three reaction-path problems with the same FAST-VQE calculator and returns the per-geometry total energies.
Plot the PES
energies = reaction_result.total_energies.values
pe_energies = pe_reaction_result.total_energies.values
even_handed_energies = even_handed_reaction_result.total_energies.values
# Plot and save to dist/
outdir = Path("dist")
plot_relative_energy_vs_coordinate(
reaction_coordinates=oh_distances,
energies_hartree=energies,
pe_energies_hartree=pe_energies,
even_handed_energies_hartree=even_handed_energies,
outdir=outdir,
)
We plot relative energies vs O–H distance for all three paths. By referencing each curve to its minimum, the shape of the PES becomes directly comparable. You should observe the Even-Handed curve is typically smoother than the geometry-wise PE curve.
Running the example
After saving the script as methanol_oh_dissociation.py, run:
python methanol_oh_dissociation.py
Results
You should see a plot:
Notes and good practice
Even-Handed PE is strongly recommended whenever a PE workflow is applied to a reaction path or any scan with large geometry changes; it stabilizes the embedded set and removes PES artifacts.
Input flexibility — instead of generating geometries in code, you can give a multi-image XYZ file (one geometry after another) and pass that file path into build_reaction_configuration(...).
References
Even-Handed localized orbital partitioning for reaction paths: https://arxiv.org/abs/1809.03004
Qrunch Projective-Embedding overview: Projection-Based Wavefunction-in-DFT Embedding