diff --git a/examples/sctep.py b/examples/sctep.py
index 404d26b..05cf229 100644
--- a/examples/sctep.py
+++ b/examples/sctep.py
@@ -1,11 +1,36 @@
 #! /usr/bin/env python
 
-from radicalpy.data import Molecule
+from radicalpy.data import Isotope, Molecule, Nucleus
+from radicalpy.experiments import steady_state_mary
+from radicalpy.simulation import Basis, LiouvilleSimulation
 
 
 def main():
-    m = Molecule(nuclei=[])
-    print(m)
+    gamma = Isotope("E").gamma_mT
+    D = -6.2 * gamma
+    E = 35 * gamma
+    triplet = Nucleus(
+        magnetogyric_ratio=gamma,
+        multiplicity=3,
+        hfc=0.0,
+        name="Triplet",
+    )
+    m = Molecule(radical=triplet)
+    sim = LiouvilleSimulation(molecules=[m, m], basis=Basis.ZEEMAN)
+    # print(sim)
+    steady_state_mary(
+        sim,
+        obs_state=State.TP_SINGLET,
+        B=Bs,
+        D=D,
+        E=E,
+        J=0,
+        kinetics=[
+            rp.kinetics.Haberkorn(krec, State.TP_SINGLET),
+            rp.kinetics.HaberkornFree(kesc),
+        ],
+    )
+    print(H)
 
 
 if __name__ == "__main__":
diff --git a/radicalpy/simulation.py b/radicalpy/simulation.py
index 3dec627..01c432f 100644
--- a/radicalpy/simulation.py
+++ b/radicalpy/simulation.py
@@ -571,22 +571,25 @@ def dipolar_hamiltonian_3d(self, dipolar_tensor: np.ndarray) -> np.ndarray:
             )
         )
 
-    def zero_field_splitting_hamiltonian(self, E, D) -> np.ndarray:
+    def zero_field_splitting_hamiltonian(self, D, E) -> np.ndarray:
         """Construct the Zero Field Splitting (ZFS) Hamiltoninan."""
-        for idx in range(self.particles):
-            spinops = []
-            coeffs = [E - D / 3, -E - D / 3, 2 * D / 3]
-            sum(
-                coeff * np.linalg.matrix_power(self.spin_operator(idx, ax), 2)
-                for idx, (coeff, ax) in itertools.product(
-                    self.particles, zip(coeffs, "xyz")
-                )
-            )
+        coeffs = [E - D / 3, -E - D / 3, 2 * D / 3]
+        result = complex(0.0)
+        other = complex(0.0)
+        for idx, p in enumerate(self.particles):
+            tmp = complex(0.0)
+            for coeff, ax in zip(coeffs, "xyz"):
+                sqr = self.spin_operator(idx, ax) @ self.spin_operator(idx, ax)
+                result += coeff * sqr
+                tmp += coeff * sqr
             Sx = self.spin_operator(idx, "x") @ self.spin_operator(idx, "x")
-            # Sx, Sy, Sz = spinops(pos, 2, spin=3)
-            # Ssquared = Sx @ Sx + Sy @ Sy + Sz @ Sz
-            # return D * (Sz @ Sz - (1 / 3) * Ssquared) + E * ((Sx @ Sx) - (Sy @ Sy))
-        return 0
+            Sy = self.spin_operator(idx, "y") @ self.spin_operator(idx, "y")
+            Sz = self.spin_operator(idx, "z") @ self.spin_operator(idx, "z")
+            Ssquared = Sx @ Sx + Sy @ Sy + Sz @ Sz
+            other = D * (Sz @ Sz - (1 / 3) * Ssquared) + E * ((Sx @ Sx) - (Sy @ Sy))
+            delta = np.linalg.norm(other - sqr)
+            print(f"{delta=}")
+        return result
 
     def total_hamiltonian(
         self,