Merge pull request #5 from molssi-seamm/dev

Significant improvements!

HISTORY.rst

@@ -1,6 +1,18 @@
 =======
 History
 =======
+2024.7.21 -- Significant improvements!
+    * Simplified error analysis to safe approach of analyzing the diffusion constants
+      over runs.
+    * Improved fitting of the curves to focus on the central linear portion. There are
+      reasonable defaults but the user can adjust as needed.
+    * Provided a combined average and error bars when both the MSF approach and Helfand
+      moments are used.
+    * Capture temperature, pressure, and cell size from the MD step, providing 1/L as a
+      result since the true diffusion constants are found by extrapolating to 1/L = 0.
+    * Provided control over the number of steps for the expensive numerical integration
+      in the Helfand moments, providing a reasonable default of 1000.
+
 2024.7.15 -- Bugfix: Significant error in Helfand Moment approach
     * Now fixed and seems to be working.
 

diffusivity_step/analysis.py

@@ -94,7 +94,7 @@ def create_helfand_moments(v, species, m=None):
 
     n, nmols, _ = v.shape
     if m is None:
-        m = min(n // 20, 1000)
+        m = min(n // 2, 5000)
 
     Ms = []
     M_errs = []
@@ -129,70 +129,7 @@ def create_helfand_moments(v, species, m=None):
     return Ms, M_errs
 
 
-def fit_helfand_moment(y, xs, sigma=None, start=0.1, end=0.9):
-    """Find the best linear fit.
-
-    Parameters
-    ----------
-    y : [float] or numpy.ndarray()
-        The Helfand moment
-
-    xs : [float]
-        The time (x) coordinate
-
-    sigma : [float] or numpy.ndarray()
-        Optional standard error of y
-
-    start : float
-        Fraction of vector to ignore at the beginning
-
-    end : float
-        The fraction of vector to end at
-
-    Returns
-    -------
-    slope : float
-        The fit slope.
-    stderr : float
-        The 95% standard error of the slope
-    xs : [float]
-        The x values (time) for the fit curve
-    ys : [float]
-        The y values for the fit curve.
-    """
-    # We know the curves curve near the origin, so ignore the first part and last
-    n = len(y)
-    i = int(n * start)
-    j = int(n * end)
-
-    if sigma is None:
-        popt, pcov, infodict, msg, ierr = curve_fit(
-            axb,
-            xs[i:j],
-            y[i:j],
-            full_output=True,
-        )
-    else:
-        popt, pcov, infodict, msg, ierr = curve_fit(
-            axb,
-            xs[i:j],
-            y[i:j],
-            full_output=True,
-            sigma=sigma[i:j],
-            absolute_sigma=True,
-        )
-    slope = float(popt[0])
-    b = float(popt[1])
-    err = float(np.sqrt(np.diag(pcov)[0]))
-
-    ys = []
-    for x in xs[i:j]:
-        ys.append(axb(x, slope, b))
-
-    return slope, err, xs[i:j], ys
-
-
-def fit_msd(y, xs, sigma=None, start=0.1, end=0.9):
+def fit_slope(y, xs, sigma=None, start=0.1, end=0.9):
     """Find the best linear fit to longest possible segment.
 
     Parameters
@@ -295,7 +232,7 @@ def add_helfand_trace(
         label, color, colora = tensor_labels[i]
         if fit is not None:
             hover = (
-                f"{species} D{label} = ({fit[i]['D_s']} ± {fit[i]['err_s']}) * "
+                f"{species} D{label} = {fit[i]['D_s']} * "
                 f"{fit[i]['scale']:.1e} m^2/s"
             )
             plot.add_trace(
@@ -466,7 +403,7 @@ def add_msd_trace(
         label, color, colora = tensor_labels[i]
         if fit is not None:
             hover = (
-                f"{species} D{label} = ({fit[i]['D_s']} ± {fit[i]['err_s']}) * "
+                f"{species} D{label} = {fit[i]['D_s']} * "
                 f"{fit[i]['scale']:.1e} m^2/s"
             )
             plot.add_trace(

diffusivity_step/data/properties.csv

@@ -8,10 +8,18 @@ Dy#MSD#{model},float,m^2/s,"The diffusion coefficient in y from MSD using {model
 Dz#MSD#{model},float,m^2/s,"The diffusion coefficient in z from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
 "Dz,stderr#MSD#{model}",float,m^2/s,"The standard error of the diffusion coefficient in z from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
 D#Helfand Moments#{model},float,m^2/s,"The total diffusion coefficient from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
-"D,stderr#Helfand Moments#{model}",float,m^2/s,"The stadnard error of the total diffusion coefficient from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
-Dx#Helfand Moments#{model},float,m^2/s,"The diffusion coefficient in x from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
-"Dx,stderr#Helfand Moments#{model}",float,m^2/s,"The standard error if the diffusion coefficient in x from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
-Dy#Helfand Moments#{model},float,m^2/s,"The diffusion coefficient in y from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
-"Dy,stderr#Helfand Moments#{model}",float,m^2/s,"The standard error if the diffusion coefficient in y from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
-Dz#Helfand Moments#{model},float,m^2/s,"The diffusion coefficient in z from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
-"Dz,stderr#Helfand Moments#{model}",float,m^2/s,"The standard error if the diffusion coefficient in z from MSD using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+"D,stderr#Helfand Moments#{model}",float,m^2/s,"The stadnard error of the total diffusion coefficient from Helfand Moments using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+Dx#Helfand Moments#{model},float,m^2/s,"The diffusion coefficient in x from Helfand Moments using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+"Dx,stderr#Helfand Moments#{model}",float,m^2/s,"The standard error of the diffusion coefficient in x from Helfand Moments using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+Dy#Helfand Moments#{model},float,m^2/s,"The diffusion coefficient in y from Helfand Moments using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+"Dy,stderr#Helfand Moments#{model}",float,m^2/s,"The standard error of the diffusion coefficient in y from Helfand Moments using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+Dz#Helfand Moments#{model},float,m^2/s,"The diffusion coefficient in z from Helfand Moments using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+"Dz,stderr#Helfand Moments#{model}",float,m^2/s,"The standard error of the diffusion coefficient in z from Helfand Moments using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+D#{model},float,m^2/s,"The total diffusion coefficient using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+"D,stderr#{model}",float,m^2/s,"The stadnard error of the total diffusion coefficient using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+Dx#{model},float,m^2/s,"The diffusion coefficient in x using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+"Dx,stderr#{model}",float,m^2/s,"The standard error of the diffusion coefficient in x using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+Dy#{model},float,m^2/s,"The diffusion coefficient in y using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+"Dy,stderr#{model}",float,m^2/s,"The standard error of the diffusion coefficient in y using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+Dz#{model},float,m^2/s,"The diffusion coefficient in z using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity
+"Dz,stderr#{model}",float,m^2/s,"The standard error of the diffusion coefficient in z using {model}",https://en.wikipedia.org/wiki/Mass_diffusivity