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