calculate_exchange_tensor based on grogu

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
+    energies = pair["energies"]
-    # dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),
+    # Initialize output arrays
-    # dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),
+    J = np.zeros((3, 3))
-    # dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),
+    D = np.zeros(3)
-
+
-    J_ii = np.array([o2[-1], o1[0], o1[-1]])  # xx, yy, zz
+    # J matrix calculations
-    J_S = -0.5 * np.array([o3[1] + o3[2], o2[1] + o2[1], o1[1] + o1[2]])  # yz, zx, xy
+    # J(1,1) = mean([DEij(2,2,2), DEij(2,2,3)])
-    D = 0.5 * np.array([o1[1] - o1[2], o2[2] - o2[1], o3[1] - o3[2]])  # x, y, z
+    J[0, 0] = np.mean([energies[1, 3], energies[2, 3]])
-    return J_ii.sum() / 3, np.concatenate([J_ii[:2] - J_ii.sum() / 3, J_S]).flatten(), D
+
+    # 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):