Tutorial: GP-Accelerated Saddle Point Search ¶

1 Finding a Transition State on the Muller-Brown Surface ¶

This tutorial walks through a complete GP-accelerated saddle point search using gpr_optim. We find the transition state between two minima on the Muller-Brown potential, a standard 2D benchmark surface.

1.1 What you will learn ¶

How to set up a GP dimer search from Python
How the GP surrogate reduces oracle (energy/force) calls
How to interpret the results

1.2 Prerequisites ¶

pip install gpr_optim ase

1.3 Step 1: Define the potential ¶

The Muller-Brown surface has three minima and two saddle points:

import numpy as np
from gpr_optim import (
    ase_potential, run_dimer, preset_parameters,
    AtomsConfiguration, Observation, PotentialWrapper,
    atoms_config_from_ase,
)

# Use the built-in Muller-Brown potential as an ASE calculator
# (In practice, replace with VASP, CP2K, or any ASE calculator)
from ase import Atoms
from ase.calculators.emt import EMT

# For this tutorial, we wrap a simple 2-atom system
atoms = Atoms("Cu2", positions=[[0, 0, 0], [2.5, 0, 0]])
atoms.calc = EMT()
pot_fn = ase_potential(atoms)

1.4 Step 2: Set up parameters ¶

# Start from the "fast" preset (fewer iterations, loose convergence)
params = preset_parameters("fast")

# Override key parameters for this system
params["dimer_sep"] = 0.01        # dimer separation
params["T_dimer"] = 0.05          # force convergence threshold
params["method_rot"] = "LBFGS_alg"
params["method_trans"] = "LBFGS_alg"
params["param_trans"] = np.array([0.1, 0.2])
params["cell_dimensions"] = np.zeros(9)

1.5 Step 3: Create initial observations ¶

# Initial midpoint observation
init_mid = Observation()
init_mid.R = np.array([[0, 0, 0, 2.5, 0, 0]], dtype=np.float64)
init_mid.E = np.array([[atoms.get_potential_energy()]], dtype=np.float64)
init_mid.G = np.array([[-atoms.get_forces().flatten()]], dtype=np.float64)

# Initial observation set (can be empty)
init_obs = Observation()
init_obs.clear()

# Initial dimer orientation
orient = np.array([[0, 0, 1, 0, 0, -1]], dtype=np.float64)
orient /= np.linalg.norm(orient)

1.6 Step 4: Set up atomic configuration ¶

atoms_config = atoms_config_from_ase(atoms)

1.7 Step 5: Run the dimer search ¶

dimer = run_dimer(params, init_obs, init_mid, orient, atoms_config, pot_fn)

print(f"Converged: {dimer.get_total_force_calls()} force calls")
print(f"Final energy: {dimer.get_final_energy():.4f}")
print(f"GPR surrogate calls: {dimer.get_total_gpr_force_calls()}")
print(f"Iterations: {dimer.get_iterations()}")

1.8 Step 6: Interpret results ¶

The key metric is the ratio of GPR surrogate calls to true force calls. A ratio > 5 means the GP saved significant computation. For expensive DFT calculations, this translates directly to wall-time savings.

1.9 Next steps ¶

Use preset_parameters("accurate") for tighter convergence
Try the T-process for robustness: see How to use the T-process
Run GP-NEB for minimum energy paths: see How to run GP-NEB