Merge branch release-2018

[alexxy/gromacs.git] / src / gromacs / tables / cubicsplinetable.cpp

diff --git a/src/gromacs/tables/cubicsplinetable.cpp b/src/gromacs/tables/cubicsplinetable.cpp
index 502b9bcb841de2806b5c52cd5d2b6810364c4473..da94daf08cd2e61cf329938a86f5569d6158b419 100644 (file)
--- a/src/gromacs/tables/cubicsplinetable.cpp
+++ b/src/gromacs/tables/cubicsplinetable.cpp
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2016,2017, by the GROMACS development team, led by
+ * Copyright (c) 2016,2017,2018, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -340,10 +340,14 @@ CubicSplineTable::CubicSplineTable(std::initializer_list<AnalyticalSplineTableIn
             ex.prependContext("Error generating cubic spline table for function '" + thisFuncInput.desc + "'");
             throw;
         }
-        // Calculate the required table spacing h. The error we make with linear interpolation
-        // of the derivative will be described by the third-derivative correction term.
-        // This means we can compute the required spacing as h = sqrt(12*tolerance*min(f'/f''')),
-        // where f'/f''' is the first and third derivative of the function, respectively.
+        // Calculate the required table spacing h. The error we make with a third order polynomial
+        // (second order for derivative) will be described by the fourth-derivative correction term.
+        //
+        // This means we can compute the required spacing as h = 0.5*cbrt(72*sqrt(3)*tolerance**min(f'/f'''')),
+        // where f'/f'''' is the first and fourth derivative of the function, respectively.
+        // Since we already have an analytical form of the derivative, we reduce the numerical
+        // errors by calculating the quotient of the function and third derivative of the
+        // input-derivative-analytical function instead.
 
         double thisMinQuotient = internal::findSmallestQuotientOfFunctionAndThirdDerivative(thisFuncInput.derivative, range_);