diff --git a/src/grogu_magn/useful.py b/src/grogu_magn/useful.py
index d3fca9a..c3db49b 100644
--- a/src/grogu_magn/useful.py
+++ b/src/grogu_magn/useful.py
@@ -216,15 +216,47 @@ def blow_up_orbindx(orb_indices):
 
 
 def calculate_exchange_tensor(pair):
-    o1, o2, o3 = pair["energies"]  # o1=x, o2=y, o3=z
-    # dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),
-    # dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),
-    # dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),
-
-    J_ii = np.array([o2[-1], o1[0], o1[-1]])  # xx, yy, zz
-    J_S = -0.5 * np.array([o3[1] + o3[2], o2[1] + o2[1], o1[1] + o1[2]])  # yz, zx, xy
-    D = 0.5 * np.array([o1[1] - o1[2], o2[2] - o2[1], o3[1] - o3[2]])  # x, y, z
-    return J_ii.sum() / 3, np.concatenate([J_ii[:2] - J_ii.sum() / 3, J_S]).flatten(), D
+    energies = pair["energies"]
+    # Initialize output arrays
+    J = np.zeros((3, 3))
+    D = np.zeros(3)
+
+    # J matrix calculations
+    # J(1,1) = mean([DEij(2,2,2), DEij(2,2,3)])
+    J[0, 0] = np.mean([energies[1, 3], energies[2, 3]])
+
+    # J(1,2) = -mean([DEij(1,2,3), DEij(2,1,3)])
+    J[0, 1] = -np.mean([energies[2, 1], energies[2, 2]])
+    J[1, 0] = J[0, 1]
+
+    # J(1,3) = -mean([DEij(1,2,2), DEij(2,1,2)])
+    J[0, 2] = -np.mean([energies[1, 1], energies[1, 2]])
+    J[2, 0] = J[0, 2]
+
+    # J(2,2) = mean([DEij(2,2,1), DEij(1,1,3)])
+    J[1, 1] = np.mean([energies[0, 3], energies[2, 0]])
+
+    # J(2,3) = -mean([DEij(1,2,1), DEij(2,1,1)])
+    J[1, 2] = -np.mean([energies[0, 1], energies[0, 2]])
+    J[2, 1] = J[1, 2]
+
+    # J(3,3) = mean([DEij(1,1,1), DEij(1,1,2)])
+    J[2, 2] = np.mean([energies[0, 0], energies[1, 0]])
+
+    # D vector calculations
+    # D(1) = mean([DEij(1,2,1), -DEij(2,1,1)])
+    D[0] = np.mean([energies[0, 1], -energies[0, 2]])
+
+    # D(2) = mean([DEij(2,1,2), -DEij(1,2,2)])
+    D[1] = np.mean([energies[1, 2], -energies[1, 1]])
+
+    # D(3) = mean([DEij(1,2,3), -DEij(2,1,3)])
+    D[2] = np.mean([energies[2, 1], -energies[2, 2]])
+
+    J_iso = np.trace(J) / 3
+    J_S = (J - J_iso * np.eye(3)).flatten()
+
+    return J_iso, J_S, D
 
 
 def print_parameters(simulation_parameters):