diff --git a/examples/iron_harmonic_fcc_bcc.py b/examples/iron_harmonic_fcc_bcc.py
index c69c5d8..e765f2c 100644
--- a/examples/iron_harmonic_fcc_bcc.py
+++ b/examples/iron_harmonic_fcc_bcc.py
@@ -4,7 +4,7 @@
 
 import psiflow
 from psiflow.data import Dataset
-from psiflow.free_energy import Integration, compute_harmonic
+from psiflow.free_energy import Integration, compute_harmonic, harmonic_free_energy
 from psiflow.geometry import Geometry
 from psiflow.hamiltonians import Harmonic, MACEHamiltonian
 from psiflow.sampling import optimize
@@ -39,8 +39,7 @@ def main():
         integration = Integration(
             harmonic,
             temperatures=[temperature],
-            delta_hamiltonian=(scaling - 1)
-            * harmonic,  # makes sampling distribution slightly broader
+            delta_hamiltonian=(scaling - 1) * harmonic,
             delta_coefficients=np.linspace(0, 1, num=4, endpoint=True),
         )
         walkers = integration.create_walkers(  # noqa: F841
@@ -51,15 +50,18 @@ def main():
         integration.compute_gradients()
 
         reduced_f = integration.along_delta(temperature=temperature).result()
-        f_harmonic = harmonic.compute_free_energy(
-            temperature=temperature, quantum=False
+        f_harmonic = harmonic_free_energy(
+            hessian,
+            temperature=temperature,
+            quantum=False,
         )
         simulated[name] = (f_harmonic.result() + reduced_f[-1]) / beta
 
         # theoretical
-        harmonic = Harmonic(minimum, scaling * hessian.result())
-        f_harmonic_scaled = harmonic.compute_free_energy(
-            temperature=temperature, quantum=False
+        f_harmonic_scaled = harmonic_free_energy(
+            scaling * hessian.result(),
+            temperature=temperature,
+            quantum=False,
         )
         theoretical[name] = f_harmonic_scaled.result() / beta
 
diff --git a/examples/water_cp2k_noise.py b/examples/water_cp2k_noise.py
index 45148b1..ef7b472 100644
--- a/examples/water_cp2k_noise.py
+++ b/examples/water_cp2k_noise.py
@@ -5,7 +5,7 @@
 import psiflow
 from psiflow.data import Dataset
 from psiflow.geometry import Geometry
-from psiflow.reference import CP2K
+from psiflow.reference import CP2K, evaluate
 
 
 def main():
@@ -42,7 +42,7 @@ def main():
     for i in range(50):
         g = copy.deepcopy(geometry)
         g.per_atom.positions += i * delta_x
-        states.append(cp2k.evaluate(g))
+        states.append(evaluate(g, cp2k))
     data = Dataset(states)
     energy = data.get("per_atom_energy").result()
     forces = data.get("forces").result()