manage units

Author: simongravelle

nmrdfrommd/NMR.py

@@ -8,6 +8,7 @@
 # Released under the GNU Public Licence, v3 or any higher version
 # SPDX-License-Identifier: GPL-3.0-or-later
 
+
 import numpy as np
 from scipy import constants as cst
 from scipy.interpolate import interp1d

nmrdfrommd/fourier.py

@@ -9,7 +9,7 @@
 # SPDX-License-Identifier: GPL-3.0-or-later
 
 import numpy as np
-from scipy import constants as cst
+from .units import PS_TO_S, HZ_TO_MHZ
 
 # This Fourier transform implementation is numerically correct
 # but scientifically fragile due to implicit unit handling and
@@ -64,10 +64,10 @@ def fourier_transform(data, normalize=False):
     if not np.allclose(np.diff(data[:, 0]), dt_ps):
         raise ValueError("Non-uniform time spacing detected.")
 
-    dt = dt_ps * cst.pico # convert to seconds
+    dt = dt_ps * PS_TO_S
 
     n = len(data)
-    freqs = np.fft.rfftfreq(n, dt) / cst.mega  # in MHz
+    freqs = np.fft.rfftfreq(n, dt) * HZ_TO_MHZ  # in MHz
     fft_norm = 'ortho' if normalize else None
     spectrum = np.fft.rfft(data[:, 1], norm=fft_norm) * dt * 2
 

nmrdfrommd/relaxation.py

@@ -13,6 +13,7 @@
 from scipy.interpolate import interp1d
 
 from .utilities import find_nearest
+from .units import ANGSTROM_TO_M
 
 
 def compute_relaxation_rates(f, J, K, isotropic):
@@ -23,11 +24,12 @@ def compute_relaxation_rates(f, J, K, isotropic):
     Parameters
     ----------
     f : np.ndarray
-        Frequency array.
+        Frequency array, in MHz.
     J : np.ndarray
-        Spectral density, shape (dim, n_freqs); dim=1 if isotropic, else 3.
+        Spectral density, shape (dim, n_freqs), in Angstrom^-6 * s;
+        dim=1 if isotropic, else 3.
     K : float
-        Dipolar relaxation prefactor.
+        Dipolar relaxation prefactor, in m^6 * s^-2 (SI).
     isotropic : bool
         Whether to use the isotropic (dim=1) or anisotropic (dim=3) formula.
 
@@ -36,21 +38,21 @@ def compute_relaxation_rates(f, J, K, isotropic):
     R1, R2 : np.ndarray
         Relaxation rates on the same frequency grid as f.
     """
-    prefactor = K / cst.angstrom ** 6
-    idx0 = find_nearest(f, 0.0)
-    J0 = J[0]
+    prefactor = K / ANGSTROM_TO_M ** 6  # Angstrom^6 * s^-2 (K rescaled m^6 -> Angstrom^6 to match J)
+    idx0 = find_nearest(f, 0.0)         # index, dimensionless
+    J0 = J[0]                           # Angstrom^-6 * s
 
     if isotropic:
-        J02 = interp1d(f, J0, fill_value="extrapolate")(2 * f)
-        R1 = prefactor * (J0 + 4 * J02) / 6
-        R2 = prefactor * (3 / 2 * J0[idx0] + 5 / 2 * J0 + J02) / 6
+        J02 = interp1d(f, J0, fill_value="extrapolate")(2 * f)            # Angstrom^-6 * s
+        R1 = prefactor * (J0 + 4 * J02) / 6                               # s^-1
+        R2 = prefactor * (3 / 2 * J0[idx0] + 5 / 2 * J0 + J02) / 6        # s^-1
     else:
-        J1 = J[1]
-        J2 = interp1d(f, J[2], fill_value="extrapolate")(2 * f)
-        R1 = prefactor * (J1 + J2)
-        R2 = prefactor * (1 / 4) * (J0[idx0] + 10 * J1 + J2)
+        J1 = J[1]                                                         # Angstrom^-6 * s
+        J2 = interp1d(f, J[2], fill_value="extrapolate")(2 * f)           # Angstrom^-6 * s
+        R1 = prefactor * (J1 + J2)                                        # s^-1
+        R2 = prefactor * (1 / 4) * (J0[idx0] + 10 * J1 + J2)              # s^-1
 
-    return R1, R2
+    return R1, R2  # s^-1, s^-1
 
 def compute_relaxation_times(f, R1, R2, target_frequency=None):
     """Compute T1 and T2 relaxation times from relaxation rates R1, R2.
@@ -62,27 +64,27 @@ def compute_relaxation_times(f, R1, R2, target_frequency=None):
     Parameters
     ----------
     f : np.ndarray
-        Frequency array.
+        Frequency array, in MHz.
     R1 : np.ndarray
-        Longitudinal relaxation rate, same shape as f.
+        Longitudinal relaxation rate, in s^-1, same shape as f.
     R2 : np.ndarray
-        Transverse relaxation rate, same shape as f.
+        Transverse relaxation rate, in s^-1, same shape as f.
     target_frequency : float, optional
-        Frequency at which to evaluate T1/T2. If None, f=0 is used.
+        Frequency at which to evaluate T1/T2, in MHz. If None, f=0 is used.
 
     Returns
     -------
     T1, T2 : float
-        Relaxation times. np.inf if the corresponding rate is ~0.
+        Relaxation times, in s. np.inf if the corresponding rate is ~0.
     """
-    target = 0.0 if target_frequency is None else target_frequency
-    idx = find_nearest(f, target)
+    target = 0.0 if target_frequency is None else target_frequency  # MHz
+    idx = find_nearest(f, target)                                   # index, dimensionless
 
-    R1_val = R1[idx]
-    R2_val = R2[idx]
+    R1_val = R1[idx]  # s^-1
+    R2_val = R2[idx]  # s^-1
 
-    eps = 1e-20
-    T1 = np.inf if np.isclose(R1_val, 0.0) else 1.0 / (R1_val + eps)
-    T2 = np.inf if np.isclose(R2_val, 0.0) else 1.0 / (R2_val + eps)
+    eps = 1e-20  # s^-1 (safety floor, sized for R1/R2 ~ 0.1-1e4 s^-1)
+    T1 = np.inf if np.isclose(R1_val, 0.0) else 1.0 / (R1_val + eps)  # s
+    T2 = np.inf if np.isclose(R2_val, 0.0) else 1.0 / (R2_val + eps)  # s
 
-    return T1, T2
+    return T1, T2  # s, s

nmrdfrommd/units.py

@@ -0,0 +1,14 @@
+"""Unit conventions and conversion factors for NMRDfromMD.
+
+Internal convention (matches MDAnalysis and common NMR usage):
+    length    : Angstrom (Å)    -- matches MDAnalysis .positions
+    time      : picosecond (ps) -- matches MDAnalysis .trajectory.dt
+    frequency : MHz             -- conventional for NMR relaxation
+Anything else (K, gamma, ...) is SI, converted to/from the above only
+at the points imported from here.
+"""
+from scipy import constants as cst
+
+PS_TO_S = cst.pico
+HZ_TO_MHZ = 1 / cst.mega
+ANGSTROM_TO_M = cst.angstrom