Merge pull request #4 from paulsaxe/main

Fixing bugs in Helfand moment fit
molssi-seamm · Jan 5, 2024 · 176f530 · 176f530
2 parents 6df3a9a + e56824e
commit 176f530
Show file tree

Hide file tree

Showing 3 changed files with 65 additions and 28 deletions.
diff --git a/HISTORY.rst b/HISTORY.rst
@@ -1,7 +1,11 @@
 =======
 History
 =======
-
+2024.1.5 -- Bugfix: Thermal conductivity
+    * If the Helfand moments fit in the thermal conductivity step failed it stopped the
+      entire job. This is correct, as well as some of the undelying causes for
+      convergence issues.
+
 2023.5.29 -- Converged with general approach for trajectory analysis
 
 2023.5.6 -- Bugfix

diff --git a/thermal_conductivity_step/analysis.py b/thermal_conductivity_step/analysis.py
@@ -1,13 +1,17 @@
 # -*- coding: utf-8 -*-
 
 """Routines to help do Green-Kubo and Helfand moments analysis."""
+import logging
+import pprint
 import warnings
 
 import numpy as np
 from scipy.integrate import cumulative_trapezoid
 from scipy.optimize import curve_fit, OptimizeWarning
 import statsmodels.tsa.stattools as stattools
 
+logger = logging.getLogger("thermal_conductivity")
+
 tensor_labels = [
     ("xx", "red", "rgba(255,0,0,0.1)"),
     ("yy", "green", "rgba(0,255,0,0.1)"),
@@ -283,7 +287,7 @@ def fit_green_kubo_integral(y, xs, sigma=None):
     return kappa, kappa_err, a, a_err, tau, tau_err, xs[1:nf], fy
 
 
-def fit_helfand_moment(y, xs, sigma=None, start=1):
+def fit_helfand_moment(y, xs, sigma=None, start=1, logger=logger):
     """Find the best linear fit to longest possible segment.
 
     Parameters
@@ -313,7 +317,7 @@ def fit_helfand_moment(y, xs, sigma=None, start=1):
     # We know the curves curve near the origin, so ignore the first part
     i = int(start / dx)
     if i > len(y):
-        i = int(1 / dx)
+        i = len(y) // 2
 
     popt, pcov, infodict, msg, ierr = curve_fit(
         axb,
@@ -323,6 +327,14 @@ def fit_helfand_moment(y, xs, sigma=None, start=1):
         sigma=sigma[i:],
         absolute_sigma=True,
     )
+    if logger.isEnabledFor(logging.DEBUG):
+        logger.debug("")
+        logger.debug(f"{popt=}")
+        logger.debug(f"{pcov=}")
+        logger.debug(pprint.pformat(infodict, compact=True))
+        logger.debug(f"{msg=}")
+        logger.debug(f"{ierr=}")
+
     slope = float(popt[0])
     b = float(popt[1])
     err = float(np.sqrt(np.diag(pcov)[0]))

diff --git a/thermal_conductivity_step/thermal_conductivity.py b/thermal_conductivity_step/thermal_conductivity.py
@@ -250,7 +250,7 @@ def description_text(self, P=None, short=False):
                         f"Error describing thermal_conductivity flowchart: {e} in "
                         f"{node}"
                     )
-                    logger.critical(
+                    self.logger.critical(
                         f"Error describing thermal_conductivity flowchart: {e} in "
                         f"{node}"
                     )
@@ -260,7 +260,7 @@ def description_text(self, P=None, short=False):
                         "Unexpected error describing thermal_conductivity flowchart: "
                         f"{sys.exc_info()[0]} in {str(node)}"
                     )
-                    logger.critical(
+                    self.logger.critical(
                         "Unexpected error describing thermal_conductivity flowchart: "
                         f"{sys.exc_info()[0]} in {str(node)}"
                     )
@@ -521,29 +521,50 @@ def analyze(self, indent="", P=None, style="full", run=None, **kwargs):
                 start = max(fit0[i]["tau"]) * 3
             else:
                 start = 1
-            slope, err, xs, ys = fit_helfand_moment(
-                self.M[i], ts, sigma=self.M_err[i], start=start
-            )
-            fit.append(
-                {
-                    "kappa": slope,
-                    "stderr": err,
-                    "xs": xs,
-                    "ys": ys,
-                }
-            )
-            v, e = fmt_err(slope, 2 * err)
-            if style == "1-line":
-                alpha = self.tensor_labels[i][0]
-                table["K" + alpha].append(v)
-                table["e" + alpha].append(e)
-            elif style == "full":
-                table["Method"].append("Helfand Moments" if i == 0 else "")
-                table["Dir"].append(self.tensor_labels[i][0])
-                table["Kappa"].append(v)
-                table["±"].append("±")
-                table["95%"].append(e)
-                table["Units"].append("W/m/K" if i == 0 else "")
+            if start < 1:
+                start = 1
+
+            if self.logger.isEnabledFor(logging.DEBUG):
+                self.logger.debug("\n\n\n**********************\n")
+                self.logger.debug(f"{i=}")
+                for v1, v2, v3 in zip(ts[0:9], self.M[i][0:9], self.M_err[i][0:9]):
+                    self.logger.debug(f"   {v1:.3f} {v2:12.4e} {v3:12.4e}")
+                self.logger.debug("...")
+                for v1, v2, v3 in zip(
+                    ts[-9:-1], self.M[i][-9:-1], self.M_err[i][-9:-1]
+                ):
+                    self.logger.debug(f"   {v1:.3f} {v2:12.4e} {v3:12.4e}")
+                self.logger.debug(f"{start=}")
+                self.logger.debug("--------")
+                self.logger.debug("")
+
+            try:
+                slope, err, xs, ys = fit_helfand_moment(
+                    self.M[i], ts, sigma=self.M_err[i], start=start, logger=self.logger
+                )
+                fit.append(
+                    {
+                        "kappa": slope,
+                        "stderr": err,
+                        "xs": xs,
+                        "ys": ys,
+                    }
+                )
+                v, e = fmt_err(slope, 2 * err)
+                if style == "1-line":
+                    alpha = self.tensor_labels[i][0]
+                    table["K" + alpha].append(v)
+                    table["e" + alpha].append(e)
+                elif style == "full":
+                    table["Method"].append("Helfand Moments" if i == 0 else "")
+                    table["Dir"].append(self.tensor_labels[i][0])
+                    table["Kappa"].append(v)
+                    table["±"].append("±")
+                    table["95%"].append(e)
+                    table["Units"].append("W/m/K" if i == 0 else "")
+            except Exception as e:
+                logger.warning("The fit of the Helfand moments failed. Continuing...")
+                logger.warning(e)
 
         self.plot_helfand_moment(self.M, ts, M_err=self.M_err, fit=fit)