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 Kvantify 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 to avoid discontinuities/cusps in the PES (see https://arxiv.org/abs/1809.03004).
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
14qc.setup_logger(qc.LOGGER_INFO)
15AtomSymbol = Literal["C", "O", "H"]
16Atom: TypeAlias = tuple[AtomSymbol, float, float, float]
17
18
19def methanol_with_stretched_bond(new_distance: float) -> list[Atom]:
20 """
21 Provide a methanol molecule with a stretched OH bond, by translating the Hydrogen along the bond direction.
22
23 Args:
24 new_distance: Target O-H distance in Å.
25
26 """
27 atoms: list[Atom] = [
28 ("C", 0.0000, 0.0000, 0.0000),
29 ("O", 1.4300, 0.0000, 0.0000), # index 1 (O)
30 ("H", 1.9300, 0.9300, 0.0000), # index 2 (OH hydrogen) <-- this one moves
31 ("H", -0.5400, 0.9400, 0.0000),
32 ("H", -0.5400, -0.4700, 0.8150),
33 ("H", -0.5400, -0.4700, -0.8150),
34 ]
35 oxygen_index = 1
36 hydrogen_index = 2
37
38 _, x_h, y_h, z_h = atoms[hydrogen_index]
39 _, x_o, y_o, z_o = atoms[oxygen_index]
40
41 vector = np.array([x_h - x_o, y_h - y_o, z_h - z_o], dtype=float)
42 vector_norm = float(np.linalg.norm(vector))
43 direction_vector = vector / vector_norm
44 new_position = np.array([x_o, y_o, z_o], dtype=float) + new_distance * direction_vector
45
46 stretched_atoms = list(atoms)
47 new_h: Atom = ("H", float(new_position[0]), float(new_position[1]), float(new_position[2]))
48 stretched_atoms[hydrogen_index] = new_h
49 return stretched_atoms
50
51
52
53
54def plot_relative_energy_vs_coordinate(
55 reaction_coordinates: NDArray[np.float64],
56 energies_hartree: list[float],
57 pe_energies_hartree: list[float],
58 even_handed_energies_hartree: list[float],
59 *,
60 outdir: Path = Path("dist"),
61 show: bool = False,
62) -> None:
63 """
64 Plot relative energy (kcal/mol) vs reaction coordinate and save PNG and CSV.
65
66 The relative energy is referenced to the minimum energy among the provided points.
67
68 Args:
69 reaction_coordinates: Reaction coordinates (Å).
70 energies_hartree: Total energies (Hartree) corresponding to q_values.
71 pe_energies_hartree: Total energies (Hartree) corresponding to q_values for Projected embedding.
72 even_handed_energies_hartree: Total energies (Hartree) corresponding to q_values for even handed.
73 outdir: Output directory where plots are written. Defaults to "dist".
74 show: Whether to display the figures interactively after saving.
75
76 """
77 outdir.mkdir(parents=True, exist_ok=True)
78
79 oh_distances = np.asarray(list(reaction_coordinates), dtype=float)
80 energies = np.asarray(list(energies_hartree), dtype=float)
81 pe_energies = np.asarray(list(pe_energies_hartree), dtype=float)
82 even_handed_energies = np.asarray(list(even_handed_energies_hartree), dtype=float)
83
84 # Compute relative energies (subtract min)
85 rel_simple = energies - np.min(energies)
86 rel_pe = pe_energies - np.min(pe_energies)
87 rel_even = even_handed_energies - np.min(even_handed_energies)
88
89 plt.figure(figsize=(6.0, 4.0))
90 plt.plot(oh_distances, rel_simple, "-o", label="Simple+Standard", color="darkgreen")
91 plt.plot(oh_distances, rel_pe, "--o", label="Simple+Projected-Embedding", color="green")
92 plt.plot(oh_distances, rel_even, ":o", label="Even-Handed", color="lightgreen")
93 plt.xlabel(r"O-H distance [Å]")
94 plt.ylabel("Energy [Hartree]")
95 plt.title("Methanol dissociation reaction PES")
96 plt.grid(visible=True, alpha=0.3)
97 plt.legend()
98 plt.tight_layout()
99
100 png_path = outdir / "methanol_rxn.png"
101 plt.savefig(png_path, dpi=200)
102 if show:
103 plt.show()
104 else:
105 plt.close()
106
107
108def main(oh_distances: NDArray[np.float64], data_path: Path | None = None) -> list[float]:
109 """
110 Scan methanol O-H dissociation where only the hydroxyl H moves.
111
112 Args:
113 oh_distances: Array of O-H distances [Å] to evaluate.
114
115 Returns:
116 List of total energies [Hartree] corresponding to the input distances.
117
118 """
119 if data_path is None:
120 directory = Path.cwd() / "data"
121 else:
122 directory = data_path
123
124 # =======================================================================
125 # Simple Reaction Path Problem Builder with standard ground state problem
126 # =======================================================================
127
128 ground_state_problem_builder = (
129 qc
130 .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
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
161 .problem_builder_creator() # Start creating a problem builder
162 .reaction_path() # Narrow: pick ground state problem
163 .simple( # Narrow: pick simple reaction path problem
164 ground_state_problem_builder # Use the ground state problem builder for each geometry
165 )
166 .create()
167 )
168
169 # ===================================================================================
170 # Simple Reaction Path Problem Builder with projective embedding ground state problem
171 # ===================================================================================
172
173 # Build Projective embedding ground state problem builder.
174 pe_ground_state_problem_builder = (
175 qc
176 .problem_builder_creator() # Start creating a problem builder
177 .ground_state() # Narrow to a ground state problem
178 .projective_embedding() # Narrow to a Projective embedding ground state problem
179 .choose_embedded_orbital_calculator()
180 .hartree_fock(
181 options=qc.options.HartreeFockCalculatorOptions(do_density_fitting=True)
182 ) # Use Hartree-Fock to construct molecular orbitals
183 .choose_repulsion_integral_builder()
184 .resolution_of_the_identity(auxiliary_basis="weigend")
185 .choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
186 .file_persister( # Perform the sub-choice selection - Here we pick file persister
187 directory=directory, # Directory where data files will be saved
188 extension=".qdk", # File extension for the data files
189 do_save=True, # Whether to save data to files
190 load_policy="fallback", # Load data if available, otherwise compute
191 overwriting_policy="rename", # Rename files if filename already exist
192 padding_width=3, # Padding width for file numbering (001, 002, ...)
193 )
194 .add_problem_modifier()
195 .to_dense_integrals()
196 .create()
197 )
198
199 # Build the reaction path problem builder
200 # with the Projective embedding ground state problem builder.
201 pe_reaction_problem_builder = (
202 qc
203 .problem_builder_creator() # Start creating a problem builder
204 .reaction_path() # Narrow: pick ground state problem
205 .simple( # Narrow: pick simple reaction path problem
206 pe_ground_state_problem_builder
207 )
208 .create()
209 )
210
211 # ===========================================================================
212 # Even handed Reaction Path Problem Builder https://arxiv.org/abs/1809.03004
213 # ===========================================================================
214
215 # Build the reaction problem builder.
216 even_handed_reaction_problem_builder = (
217 qc
218 .problem_builder_creator() # Start creating a problem builder
219 .reaction_path() # Narrow: pick ground state problem
220 .even_handed() # Narrow: pick even-handed reaction path problem
221 .choose_embedded_orbital_calculator()
222 .hartree_fock(
223 options=qc.options.HartreeFockCalculatorOptions(do_density_fitting=True)
224 ) # Use Hartree-Fock to construct molecular orbitals
225 .choose_repulsion_integral_builder()
226 .resolution_of_the_identity(auxiliary_basis="weigend")
227 .choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
228 .file_persister( # Perform the sub-choice selection - Here we pick file persister
229 directory=directory, # Directory where data files will be saved
230 extension=".qdk", # File extension for the data files
231 do_save=True, # Whether to save data to files
232 load_policy="fallback", # Load data if available, otherwise compute
233 overwriting_policy="rename", # Rename files if filename already exist
234 padding_width=3, # Padding width for file numbering (001, 002, ...)
235 )
236 .add_problem_modifier()
237 .to_dense_integrals()
238 .create()
239 )
240
241 # Build the FAST-VQE calculator
242 fast_vqe_calculator = (
243 qc
244 .calculator_creator() # Start creating a calculator
245 .vqe() # Narrow: pick Variational quantum eigensolver (VQE)
246 .iterative() # Narrow: pick the iterative VQE
247 .standard() # Narrow: pick the standard ansatz VQE (FAST-VQE)
248 .with_options(options=qc.options.IterativeVqeOptions(max_iterations=200))
249 .choose_stopping_criterion()
250 .patience(patience=2, threshold=1e-3)
251 .choose_data_persister_manager() # Start sub-choice: Choose the data persister manager
252 .file_persister( # Perform the sub-choice selection - Here we pick file persister
253 directory=directory, # Directory where data files will be saved
254 extension=".qdk", # File extension for the data files
255 do_save=True, # Whether to save data to files
256 load_policy="fallback", # Load data if available, otherwise compute
257 overwriting_policy="rename", # Rename files if filename already exist
258 padding_width=3, # Padding width for file numbering (001, 002, ...)
259 )
260 .create() # Create the calculator instance
261 )
262
263 # Build the list of molecules that make up the reaction path.
264 reaction = []
265 for distance in oh_distances:
266 molecule = methanol_with_stretched_bond(distance)
267 reaction.append(molecule)
268
269 reaction_configuration = qc.build_reaction_configuration(
270 reaction=reaction,
271 basis_set="sto3g",
272 spin_difference=0,
273 charge=0,
274 embedded_atoms=[1, 2],
275 aux_basis_set="weigend",
276 )
277
278 reaction_path_problem = reaction_problem_builder.build_unrestricted(reaction_configuration)
279 reaction_result = fast_vqe_calculator.calculate(reaction_path_problem)
280
281 pe_reaction_path_problem = pe_reaction_problem_builder.build_unrestricted(reaction_configuration)
282 pe_reaction_result = fast_vqe_calculator.calculate(pe_reaction_path_problem)
283
284 even_handed_reaction_problem = even_handed_reaction_problem_builder.build_unrestricted(reaction_configuration)
285 even_handed_reaction_result = fast_vqe_calculator.calculate(even_handed_reaction_problem)
286
287 energies = reaction_result.total_energies.values
288 pe_energies = pe_reaction_result.total_energies.values
289 even_handed_energies = even_handed_reaction_result.total_energies.values
290
291 # Plot and save to dist/
292 outdir = Path("dist")
293 plot_relative_energy_vs_coordinate(
294 reaction_coordinates=oh_distances,
295 energies_hartree=energies,
296 pe_energies_hartree=pe_energies,
297 even_handed_energies_hartree=even_handed_energies,
298 outdir=outdir,
299 )
300
301
302 return reaction_result.total_energies.values
303
304
305if __name__ == "__main__":
306 # OH distances in Å
307 distances = np.linspace(1.0, 3.4, 6, dtype=np.float64)
308 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
The public and stable Kvantify Qrunch API is loaded as
import qrunch as qc. For most users this is the only import needed, and the one guaranteed to remain stable across minor 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 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 PE — 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
Kvantify Qrunch Projective-Embedding overview: Projection-Based Wavefunction-in-DFT Embedding