Update semiclassical_std

radicalpy/data.py

@@ -582,9 +582,6 @@ def semiclassical_std(self) -> float:
 
         Calculate the standard deviation :math:`\sigma` where
 
-        .. todo::
-           Do the math properly.
-
         .. math::
            \sigma = \sqrt{\frac{2}{\tau^2}}
 
@@ -599,14 +596,16 @@ def semiclassical_std(self) -> float:
         Examples:
             >>> m = Molecule.fromdb("flavin_anion", nuclei=["N14"])
             >>> m.semiclassical_std
-            0.0010410773656027476
+            7.663920853309001e-07
 
+        .. todo::
+           reference
         """
-        tmp = sum(
+        tau = 6 / sum(
             n.spin_quantum_number * (n.spin_quantum_number + 1) * n.hfc.isotropic**2
             for n in self.nuclei
         )
-        return tmp**2 / 6
+        return np.sqrt(2) / tau
 
 
 class Triplet(Molecule):

utils/wip_experiments.py

@@ -113,16 +113,15 @@ def algebra():
     from sympy.abc import T, h, l, n, sigma, tau, x
 
     pexpr = (tau**2 / (4 * pi)) ** (3 / 2) * exp(-1 / 4 * x**2 * tau**2)
-    pexpr = pexpr.subs(tau**2, 1 / T).subs(x, h)
+    pexpr = pexpr.subs(x, h)
     print(pexpr)
 
     bexpr = (3 / (2 * pi * n * l**2)) ** (3 / 2) * exp(-(3 * h**2) / (2 * n * l**2))
-    bexpr = bexpr.subs(n * l**2, T * 6)
+    bexpr = bexpr.subs(n * l**2, 6 * tau**-2)
     print(bexpr)
+
     pb_delta = (pexpr - bexpr).simplify()
     print(f"{pb_delta=}")
-    # sigma = 6T = 6 tau^-2
-    # 1/T = tau^2
 
     nexpr = 1 / ((sigma**2 * 2 * pi) ** (1 / 2)) * exp(-1 / 2 * x**2 / sigma**2)
     nexpr = nexpr.subs(sigma**2, 2 * T).subs(x, h)