From 062d4b608185292fcb22b6fb2d682d6f6e720ee2 Mon Sep 17 00:00:00 2001 From: David van der Spoel Date: Thu, 8 Jan 2015 11:16:27 +0100 Subject: [PATCH] Updated exponential fitting to make it robust. The exponential fitting code was not robust under different data. This is improved by ensuring that e.g. time constants are always positive. For multi-exponential fits the time constants are guaranteed to be in increasing order. Rewrote test code. Added test for error estimation that reproduces 5.0 behavior. Added new manual section "Curve fitting in GROMACS" (8.6). General cleanups of the code were done. Change-Id: Ib72fccf7f85742afeeb3fc0fd6fbd44c1c47795a --- docs/manual/analyse.tex | 111 +- docs/manual/monster.bib | 50 + src/external/lmfit/lmcurve.c | 20 +- src/external/lmfit/lmcurve.h | 2 +- src/gromacs/correlationfunctions/expfit.c | 698 ------------- src/gromacs/correlationfunctions/expfit.cpp | 986 ++++++++++++++++++ src/gromacs/correlationfunctions/expfit.h | 27 +- .../correlationfunctions/tests/expfit.cpp | 145 ++- .../tests/refdata/AutocorrTest_EacCos.xml | 42 +- .../tests/refdata/AutocorrTest_EacNormal.xml | 146 +-- .../tests/refdata/AutocorrTest_EacP0.xml | 18 +- .../tests/refdata/AutocorrTest_EacP1.xml | 10 +- .../tests/refdata/AutocorrTest_EacP2.xml | 220 ++-- .../tests/refdata/AutocorrTest_EacP4.xml | 18 +- .../tests/refdata/AutocorrTest_EacVector.xml | 18 +- .../tests/refdata/ExpfitTest_EffnERF.xml | 8 +- .../tests/refdata/ExpfitTest_EffnERREST.xml | 6 +- .../tests/refdata/ExpfitTest_EffnEXP1.xml | 2 +- .../tests/refdata/ExpfitTest_EffnEXP2.xml | 4 +- .../tests/refdata/ExpfitTest_EffnEXP5.xml | 10 +- .../tests/refdata/ExpfitTest_EffnEXP7.xml | 14 +- .../tests/refdata/ExpfitTest_EffnEXP9.xml | 18 +- ...EffnEXP3.xml => ExpfitTest_EffnEXPEXP.xml} | 6 +- .../tests/refdata/ExpfitTest_EffnPRES.xml | 12 +- .../tests/refdata/ExpfitTest_EffnVAC.xml | 4 +- .../correlationfunctions/tests/testERF.xvg | 401 +++++++ .../correlationfunctions/tests/testERREST.xvg | 501 +++++++++ .../correlationfunctions/tests/testEXP.xvg | 501 --------- .../correlationfunctions/tests/testINVEXP.xvg | 0 .../tests/testINVEXP79.xvg | 501 +++++++++ src/gromacs/gmxana/gmx_analyze.c | 127 ++- src/gromacs/gmxana/gmx_densorder.cpp | 10 +- src/gromacs/gmxana/gmx_dielectric.c | 2 +- src/gromacs/gmxana/gmx_hbond.c | 2 +- src/gromacs/gmxana/gmx_tcaf.c | 4 +- 35 files changed, 3031 insertions(+), 1613 deletions(-) delete mode 100644 src/gromacs/correlationfunctions/expfit.c create mode 100644 src/gromacs/correlationfunctions/expfit.cpp mode change 100755 => 100644 src/gromacs/correlationfunctions/tests/expfit.cpp rename src/gromacs/correlationfunctions/tests/refdata/{ExpfitTest_EffnEXP3.xml => ExpfitTest_EffnEXPEXP.xml} (63%) create mode 100644 src/gromacs/correlationfunctions/tests/testERF.xvg create mode 100644 src/gromacs/correlationfunctions/tests/testERREST.xvg delete mode 100644 src/gromacs/correlationfunctions/tests/testEXP.xvg mode change 100755 => 100644 src/gromacs/correlationfunctions/tests/testINVEXP.xvg create mode 100644 src/gromacs/correlationfunctions/tests/testINVEXP79.xvg diff --git a/docs/manual/analyse.tex b/docs/manual/analyse.tex index ff2a382f06..d235e9c454 100644 --- a/docs/manual/analyse.tex +++ b/docs/manual/analyse.tex @@ -1,7 +1,7 @@ % % This file is part of the GROMACS molecular simulation package. % -% Copyright (c) 2013,2014, by the GROMACS development team, led by +% Copyright (c) 2013,2014,2015, 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. @@ -479,6 +479,115 @@ but this is not very accurate~\cite{Hess2002a}, and actually the values do not converge. %} % Brace matches ifthenelse test for gmxlite +\section{Curve fitting in \gromacs} +\subsection{Sum of exponential functions} +Sometimes it is useful to fit a curve to an analytical function, for +example in the case of autocorrelation functions with noisy +tails. {\gromacs} is not a general purpose curve-fitting tool however +and therefore {\gromacs} only supports a limited number of +functions. +Table~\ref{tab:fitfn} lists the available options with the +corresponding command-line options. The underlying routines for +fitting use the Levenberg-Marquardt algorithm as implemented in the +{\tt lmfit} package~\cite{lmfit} (a bare-bones version of which is +included in {\gromacs} in which an option for error-weighted fitting +was implemented). +\begin{table}[ht] +\centering +\caption{Overview of fitting functions supported in (most) analysis tools + that compute autocorrelation functions. The ``Note'' column describes + properties of the output parameters.} +\label{tab:fitfn} +\begin{tabular}{lcc} +\hline +Command & Functional form $f(t)$& Note\\ +line option & & \\ +\hline +exp & $e^{-t/{a_0}}$ &\\ +aexp & $a_1e^{-t/{a_0}}$ &\\ +exp_exp & $a_1e^{-t/{a_0}}+(1-a_1)e^{-t/{a_2}}$ & $a_2\ge a_0\ge 0$\\ +exp5 & $a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_4$ &$a_2\ge a_0\ge 0$\\ +exp7 & +$a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_5e^{-t/{a_4}}+a_6$&$a_4\ge a_2\ge +a_0 \ge0$\\ +exp9 & +$a_1e^{-t/{a_0}}+a_3e^{-t/{a_2}}+a_5e^{-t/{a_4}}+a_7e^{-t/{a_6}}+a_8$&$a_6\ge +a_4\ge a_2\ge a_0\ge 0$\\ +\hline +\end{tabular} +\end{table} + +\subsection{Error estimation} +Under the hood {\gromacs} implements some more fitting functions, +namely a function to estimate the error in time-correlated data due to Hess~\cite{Hess2002a}: +\beq +\varepsilon^2(t) = +\alpha\tau_1\left(1+\frac{\tau_1}{t}\left(e^{-t/\tau_1}-1\right)\right) + + (1-\alpha)\tau_2\left(1+\frac{\tau_2}{t}\left(e^{-t/\tau_2}-1\right)\right) +\eeq +where $\tau_1$ and +$\tau_2$ are time constants (with $\tau_2 \ge \tau_1$) and $\alpha$ +usually is close to 1 (in the fitting procedure it is enforced that +$0\leq\alpha\leq 1$). This is used in {\tt gmx analyze} for error +estimation using +\beq +\lim_{t\rightarrow\infty}\varepsilon(t) = \sigma\sqrt{\frac{2(\alpha\tau_1+(1-\alpha)\tau_2)}{T}} +\eeq +where $\sigma$ is the standard deviation of the data set and $T$ is +the total simulation time~\cite{Hess2002a}. + +\subsection{Interphase boundary demarcation} +In order to determine the position and width of an interface, +Steen-S{\ae}thre {\em et al.} fitted a density profile to the +following function +\beq +f(x) ~=~ \frac{a_0+a_1}{2} - \frac{a_0-a_1}{2}{\rm + erf}\left(\frac{x-a_2}{a_3^2}\right) +\eeq +where $a_0$ and $a_1$ are densities of different phases, $x$ is the +coordinate normal to the interface, $a_2$ is the position of the +interface and $a_3$ is the width of the +interface~\cite{Steen-Saethre2014a}. +This is implemented in {\tt gmx densorder}. + +\subsection{Transverse current autocorrelation function} +In order to establish the transverse current autocorrelation function +(useful for computing viscosity~\cite{Palmer1994a}) +the following function is fitted: +\beq +f(x) ~=~ e^{-\nu}\left({\rm cosh}(\omega\nu)+\frac{{\rm + sinh}(\omega\nu)}{\omega}\right) +\eeq +with $\nu = x/(2a_0)$ and $\omega = \sqrt{1-a_1}$. +This is implemented in {\tt gmx tcaf}. + +\subsection{Viscosity estimation from pressure autocorrelation + function} +The viscosity is a notoriously difficult property to extract from +simulations~\cite{Hess2002a,Wensink2003a}. It is {\em in principle} +possible to determine it by integrating the pressure autocorrelation +function~\cite{PSmith93c}, however this is often hampered by the noisy +tail of the ACF. A workaround to this is fitting the ACF to the +following function~\cite{Guo2002b}: +\beq +f(t)/f(0) = (1-C) {\rm cos}(\omega t) e^{-(t/\tau_f)^{\beta_f}} + C +e^{-(t/\tau_s)^{\beta_s}} +\eeq +where $\omega$ is the frequency of rapid pressure oscillations (mainly +due to bonded forces in molecular simulations), $\tau_f$ and $\beta_f$ +are the time constant and exponent of fast relaxation in a +stretched-exponential approximation, $\tau_s$ and $\beta_s$ are constants +for slow relaxation and $C$ is the pre-factor that determines the +weight between fast and slow relaxation. After a fit, the integral of +the function $f(t)$ is used to compute the viscosity: +\beq +\eta = \frac{V}{k_B T}\int_0^{\infty} f(t) dt +\eeq +This equation has been +applied to computing the bulk and shear viscosity using different +elements from the pressure tensor~\cite{Fanourgakis2012a}. +This is implemented in {\tt gmx viscosity}. + \section{Mean Square Displacement} \label{sec:msd} {\tt gmx msd}\\ diff --git a/docs/manual/monster.bib b/docs/manual/monster.bib index 0774ae539c..a7ba8b5b9d 100644 --- a/docs/manual/monster.bib +++ b/docs/manual/monster.bib @@ -8691,3 +8691,53 @@ pages = {105--120}, keywords = {CHARMM, force field, molecular dynamics, parametrization, DNA, RNA}, year = {2000}, } + +@Misc{lmfit, + author = {Joachim Wuttke}, + title = {lmfit}, + howpublished = {http://apps.jcns.fz-juelich.de/doku/sc/lmfit}, + year = 2013} + +@Article{Steen-Saethre2014a, + author = {B. Steen-S{\ae}thre and A. C. Hoffmann and D. van der Spoel}, + title = {Order Parameters and Algorithmic Approaches for Detection and Demarcation of Interfaces in Hydrate−Fluid and Ice−Fluid Systems}, + journal = {J. Chem. Theor. Comput.}, + year = 2014, + volume = 10, + pages = {5606-5615}} + +@Article{ Wensink2003a, + author = "E. J. W. Wensink and A. C. Hoffmann and P. J. {van + Maaren} and D. {van der Spoel}", + title = "Dynamic properties of water/alcohol mixtures studied + by computer simulation", + journal = BTjcp, + year = 2003, + volume = 119, + pages = "7308-7317" +} + +@Article{Guo2002b, + author = {Guang-Jun Guo and Yi-Gang Zhang and Keith Refson and Ya-Juan Zhao}, + title = {Viscosity and stress autocorrelation function in supercooled water: a molecular dynamics study}, + journal = {Mol. Phys.}, + year = 2002, + volume = 100, + pages = {2617-2627}} + +@Article{Fanourgakis2012a, + author = {George S. Fanourgakis and J. S. Medina and R. Prosmiti}, + title = {Determining the Bulk Viscosity of Rigid Water Models}, + journal = {J. Phys. Chem. A}, + year = 2012, + volume = 116, + pages = {2564-2570}} + +@Article{Palmer1994a, + author = {B. J. Palmer}, + title = {Transverse-current autocorrelation-function calculations of the shear viscosity for molecular liquids.}, + journal = {Phys. Rev. E}, + year = 1994, + volume = 49, + pages = {359-366}} + diff --git a/src/external/lmfit/lmcurve.c b/src/external/lmfit/lmcurve.c index 838f884bb6..4057fbff83 100644 --- a/src/external/lmfit/lmcurve.c +++ b/src/external/lmfit/lmcurve.c @@ -18,6 +18,7 @@ typedef struct { const double *t; const double *y; + const double *dy; double (*f)(double t, const double *par); } lmcurve_data_struct; @@ -25,27 +26,28 @@ typedef struct { void lmcurve_evaluate( const double *par, int m_dat, const void *data, double *fvec, gmx_unused int *info ) { - int i; + int i; + double fy; + lmcurve_data_struct *d = (lmcurve_data_struct*) data; for (i = 0; i < m_dat; i++) { - fvec[i] = - ((lmcurve_data_struct*)data)->y[i] - - ((lmcurve_data_struct*)data)->f( - ((lmcurve_data_struct*)data)->t[i], par ); + fy = d->f(d->t[i], par ); + fvec[i] = (d->y[i] - fy)/d->dy[i]; } } void lmcurve( int n_par, double *par, int m_dat, - const double *t, const double *y, + const double *t, const double *y, const double *dy, double (*f)( double t, const double *par ), const lm_control_struct *control, lm_status_struct *status ) { lmcurve_data_struct data; - data.t = t; - data.y = y; - data.f = f; + data.t = t; + data.y = y; + data.dy = dy; + data.f = f; lmmin( n_par, par, m_dat, (const void*) &data, lmcurve_evaluate, control, status ); diff --git a/src/external/lmfit/lmcurve.h b/src/external/lmfit/lmcurve.h index d74114caaf..756805e6e0 100644 --- a/src/external/lmfit/lmcurve.h +++ b/src/external/lmfit/lmcurve.h @@ -29,7 +29,7 @@ __BEGIN_DECLS void lmcurve( int n_par, double *par, int m_dat, - const double *t, const double *y, + const double *t, const double *y, const double *dy, double (*f)( double t, const double *par ), const lm_control_struct *control, lm_status_struct *status ); diff --git a/src/gromacs/correlationfunctions/expfit.c b/src/gromacs/correlationfunctions/expfit.c deleted file mode 100644 index 490e0be4f8..0000000000 --- a/src/gromacs/correlationfunctions/expfit.c +++ /dev/null @@ -1,698 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 1991-2000, University of Groningen, The Netherlands. - * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014, 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. - * - * GROMACS is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public License - * as published by the Free Software Foundation; either version 2.1 - * of the License, or (at your option) any later version. - * - * GROMACS is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with GROMACS; if not, see - * http://www.gnu.org/licenses, or write to the Free Software Foundation, - * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. - * - * If you want to redistribute modifications to GROMACS, please - * consider that scientific software is very special. Version - * control is crucial - bugs must be traceable. We will be happy to - * consider code for inclusion in the official distribution, but - * derived work must not be called official GROMACS. Details are found - * in the README & COPYING files - if they are missing, get the - * official version at http://www.gromacs.org. - * - * To help us fund GROMACS development, we humbly ask that you cite - * the research papers on the package. Check out http://www.gromacs.org. - */ -/*! \internal - * \file - * \brief - * Implements routine for fitting a data set to a curve - * - * \author David van der Spoel - * \ingroup module_correlationfunctions - */ -#include "gmxpre.h" - -#include "expfit.h" - -#include -#include - -#include "external/lmfit/lmcurve.h" - -#include "gromacs/correlationfunctions/integrate.h" -#include "gromacs/fileio/xvgr.h" -#include "gromacs/legacyheaders/macros.h" -#include "gromacs/math/vec.h" -#include "gromacs/utility/fatalerror.h" -#include "gromacs/utility/futil.h" -#include "gromacs/utility/real.h" -#include "gromacs/utility/smalloc.h" - -/*! \brief Number of parameters for each fitting function */ -static const int nfp_ffn[effnNR] = { 0, 1, 2, 3, 2, 5, 7, 9, 4, 3, 6}; - -const char *s_ffn[effnNR+2] = { - NULL, "none", "exp", "aexp", "exp_exp", "vac", - "exp5", "exp7", "exp9", "erffit", "effnPRES", NULL, NULL -}; -/* We don't allow errest as a choice on the command line */ - -/*! \brief Long description for each fitting function type */ -static const char *longs_ffn[effnNR] = { - "no fit", - "y = exp(-x/a1)", - "y = a2 exp(-x/a1)", - "y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)", - "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a1), w = sqrt(1 - a2)", - "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5", - "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7", - "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7 exp(-x/a8) + a9", - "y = 1/2*(a1+a2) - 1/2*(a1-a2)*erf( (x-a3) /a4^2)", - "y = a2 *2*a1*((exp(-x/a1)-1)*(a1/x)+1)+(1-a2)*2*a3*((exp(-x/a3)-1)*(a3/x)+1)", - "y = (1-a1) * exp(-(x/a3)^a4)*cos(x*a2) + a1*exp(-(x/a5)^a6)" -}; - -int effnNparams(int effn) -{ - if ((0 <= effn) && (effn < effnNR)) - { - return nfp_ffn[effn]; - } - else - { - return -1; - } -} - -const char *effnDescription(int effn) -{ - if ((0 <= effn) && (effn < effnNR)) - { - return longs_ffn[effn]; - } - else - { - return NULL; - } -} - -int sffn2effn(const char **sffn) -{ - int eFitFn, i; - - eFitFn = 0; - for (i = 0; i < effnNR; i++) - { - if (sffn[i+1] && strcmp(sffn[0], sffn[i+1]) == 0) - { - eFitFn = i; - } - } - - return eFitFn; -} - -/*! \brief Compute exponential function A exp(-x/tau) */ -static double myexp(double x, double A, double tau) -{ - if ((A == 0) || (tau == 0)) - { - return 0; - } - return A*exp(-x/tau); -} - -/*! \brief Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2) */ -static double lmc_erffit (double x, const double *a) -{ - double erfarg; - - erfarg = (x-a[2])/(a[3]*a[3]); - return (a[0]+a[1])/2-(a[0]-a[1])*gmx_erf(erfarg); -} - -/*! \brief Compute y = exp(-x/a0) */ -static double lmc_exp_one_parm(double x, const double *a) -{ - return exp(-x/a[0]); -} - -/*! \brief Compute y = a1 exp(-x/a0) */ -static double lmc_exp_two_parm(double x, const double *a) -{ - return a[1]*exp(-x/a[0]); -} - -/*! \brief Compute y = a1 exp(-x/a0) + (1-a1) exp(-x/a2) */ -static double lmc_exp_3_parm(double x, const double *a) -{ - double e1, e2; - - e1 = exp(-x/a[0]); - e2 = exp(-x/a[2]); - return a[1]*e1 + (1-a[1])*e2; -} - -/*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 */ -static double lmc_exp_5_parm(double x, const double *a) -{ - double e1, e2; - - e1 = exp(-x/a[1]); - e2 = exp(-x/a[3]); - return a[0]*e1 + a[2]*e2 + a[4]; -} - -/*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 exp(-x/a5) + a6 */ -static double lmc_exp_7_parm(double x, const double *a) -{ - double e1, e2, e3; - - e1 = exp(-x/a[1]); - e2 = exp(-x/a[3]); - e3 = exp(-x/a[5]); - return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6]; -} - -/*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 exp(-x/a5) + a6 exp(-x/a7) + a8 */ -static double lmc_exp_9_parm(double x, const double *a) -{ - double e1, e2, e3, e4; - - e1 = exp(-x/a[1]); - e2 = exp(-x/a[3]); - e3 = exp(-x/a[5]); - e4 = exp(-x/a[7]); - return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6]*e4 + a[8]; -} - -/*! \brief Compute y = (1-a1) * exp(-(x/a3)^a4)*cos(x*a2) + a1*exp(-(x/a5)^a6) */ -static double effnPRES_fun(double x, const double *a) -{ - double term1, term2, term3; - - if (a[0] != 0) - { - term3 = a[0] * exp(-pow((x/a[4]), a[5])); - } - else - { - term3 = 0; - } - - term1 = 1-a[0]; - if (term1 != 0) - { - term2 = exp(-pow((x/a[2]), a[3])) * cos(x*a[1]); - } - else - { - term2 = 0; - } - - return term1*term2 + term3; -} - -/*! \brief Compute vac function */ -static double lmc_vac_2_parm(double x, const double *a) -{ - /* Fit to function - * - * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v) - * - * = exp(-v) (cosh(wv) + 1/w sinh(wv)) - * - * v = x/(2 a1) - * w = sqrt(1 - a2) - * - * For tranverse current autocorrelation functions: - * a1 = tau - * a2 = 4 tau (eta/rho) k^2 - * - */ - - double y, v, det, omega, omega2, em, ec, es; - - v = x/(2*a[0]); - det = 1 - a[1]; - em = exp(-v); - if (det != 0) - { - omega2 = fabs(det); - omega = sqrt(omega2); - if (det > 0) - { - ec = em*0.5*(exp(omega*v)+exp(-omega*v)); - es = em*0.5*(exp(omega*v)-exp(-omega*v))/omega; - } - else - { - ec = em*cos(omega*v); - es = em*sin(omega*v)/omega; - } - y = ec + es; - } - else - { - y = (1+v)*em; - } - return y; -} - -/*! \brief Compute error estimate */ -static double lmc_errest_3_parm(double x, const double *a) -{ - /* - e1 = (exp(-x/a1) - 1) - e2 = (exp(-x/a3) - 1) - v1= 2*a1 * (e1*a1/x + 1) - v2 = 2*a3 * (e2*a3/x + 1) - fun = a2*v1 + (1 - a2) * v2 - */ - double e1, e2, v1, v2; - - if (a[0] != 0) - { - e1 = gmx_expm1(-x/a[0]); - } - else - { - e1 = 0; - } - if (a[2] != 0) - { - e2 = gmx_expm1(-x/a[2]); - } - else - { - e2 = 0; - } - - if (x > 0) - { - v1 = 2*a[0]*(e1*a[0]/x + 1); - v2 = 2*a[2]*(e2*a[2]/x + 1); - return a[1]*v1 + (1-a[1])*v2; - } - else - { - return 0; - } -} - -/*! \brief function type for passing to fitting routine */ -typedef double (*t_lmcurve)(double x, const double *a); - -/*! \brief array of fitting functions corresponding to the pre-defined types */ -t_lmcurve lmcurves[effnNR] = { - lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm, - lmc_exp_3_parm, lmc_vac_2_parm, - lmc_exp_5_parm, lmc_exp_7_parm, - lmc_exp_9_parm, lmc_erffit, lmc_errest_3_parm, effnPRES_fun -}; - -double fit_function(int eFitFn, double *parm, double x) -{ - // TODO: check that eFitFn is within bounds - return lmcurves[eFitFn](x, parm); -} - -/*! \brief lmfit_exp supports up to 3 parameter fitting of exponential functions */ -static gmx_bool lmfit_exp(int nfit, double x[], double y[], - real ftol, int maxiter, - double parm[], gmx_bool bVerbose, - int eFitFn) -{ - double chisq, ochisq; - gmx_bool bCont; - int i, j; - lm_control_struct *control; - lm_status_struct *status; - - if ((eFitFn < 0) || (eFitFn >= effnNR)) - { - gmx_fatal(FARGS, "fitfn = %d, should be in 0..%d (%s,%d)", - effnNR-1, eFitFn, __FILE__, __LINE__); - } - - snew(control, 1); - control->ftol = ftol; - control->xtol = ftol; - control->gtol = ftol; - control->epsilon = 0.1; - control->stepbound = 100; - control->patience = maxiter; - control->scale_diag = 1; - control->msgfile = stdout; - control->verbosity = (bVerbose ? 1 : 0); - control->n_maxpri = nfp_ffn[eFitFn]; - control->m_maxpri = 0; - snew(status, 1); - /* Initial params */ - chisq = 1e12; - j = 0; - if (bVerbose) - { - fprintf(stderr, "%4s %10s Parameters\n", - "Step", "chi^2"); - } - do - { - ochisq = chisq; - - lmcurve(nfp_ffn[eFitFn], parm, nfit, x, y, - lmcurves[eFitFn], control, status); - chisq = sqr(status->fnorm); - if (bVerbose) - { - int mmm; - fprintf(stderr, "%4d %10.5e", j, chisq); - for (mmm = 0; (mmm < nfp_ffn[eFitFn]); mmm++) - { - fprintf(stderr, " %10.5e", parm[mmm]); - } - fprintf(stderr, "\n"); - } - j++; - bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) || - ((ochisq == chisq))); - } - while (bCont && (j < maxiter)); - if (bVerbose) - { - fprintf(stderr, "\n"); - } - - sfree(control); - sfree(status); - - return TRUE; -} - -/*! \brief See description in header file. */ -real do_lmfit(int ndata, real c1[], real sig[], real dt, real x0[], - real begintimefit, real endtimefit, const output_env_t oenv, - gmx_bool bVerbose, int eFitFn, double fitparms[], int fix) -{ - FILE *fp; - char buf[32]; - - int i, j, nparm, nfitpnts; - double integral, ttt; - double *parm, *dparm; - double *x, *y, *dy; -#ifdef GMX_DOUBLE - real ftol = 1e-16; -#else - real ftol = 1e-8; -#endif - int maxiter = 50; - - if (0 != fix) - { - fprintf(stderr, "Using fixed parameters in curve fitting is not working anymore\n"); - } - nparm = nfp_ffn[eFitFn]; - if (debug) - { - fprintf(debug, "There are %d points to fit %d vars!\n", ndata, nparm); - fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n", - eFitFn, begintimefit, endtimefit, dt); - } - - snew(x, ndata); - snew(y, ndata); - snew(dy, ndata); - - j = 0; - for (i = 0; i < ndata; i++) - { - ttt = x0 ? x0[i] : dt*i; - if (ttt >= begintimefit && ttt <= endtimefit) - { - x[j] = ttt; - y[j] = c1[i]; - - /* mrqmin does not like sig to be zero */ - if (sig[i] < 1.0e-7) - { - dy[j] = 1.0e-7; - } - else - { - dy[j] = sig[i]; - } - if (debug) - { - fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy= %g\n", - j, i, x[j], y[j], dy[j]); - } - j++; - } - } - nfitpnts = j; - integral = 0; - if (nfitpnts < nparm) - { - fprintf(stderr, "Not enough data points for fitting!\n"); - } - else - { - snew(parm, nparm); - snew(dparm, nparm); - if (fitparms) - { - for (i = 0; (i < nparm); i++) - { - parm[i] = fitparms[i]; - } - } - - if (!lmfit_exp(nfitpnts, x, y, ftol, maxiter, parm, bVerbose, eFitFn)) - { - fprintf(stderr, "Fit failed!\n"); - } - else if (nparm <= 3) - { - /* Compute the integral from begintimefit */ - if (nparm == 3) - { - integral = (parm[0]*myexp(begintimefit, parm[1], parm[0]) + - parm[2]*myexp(begintimefit, 1-parm[1], parm[2])); - } - else if (nparm == 2) - { - integral = parm[0]*myexp(begintimefit, parm[1], parm[0]); - } - else if (nparm == 1) - { - integral = parm[0]*myexp(begintimefit, 1, parm[0]); - } - else - { - gmx_fatal(FARGS, "nparm = %d in file %s, line %d", - nparm, __FILE__, __LINE__); - } - - /* Generate THE output */ - if (bVerbose) - { - fprintf(stderr, "FIT: # points used in fit is: %d\n", nfitpnts); - fprintf(stderr, "FIT: %21s%21s%21s\n", - "parm0 ", "parm1 (ps) ", "parm2 (ps) "); - fprintf(stderr, "FIT: ------------------------------------------------------------\n"); - fprintf(stderr, "FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n", - parm[0], dparm[0], parm[1], dparm[1], parm[2], dparm[2]); - fprintf(stderr, "FIT: Integral (calc with fitted function) from %g ps to inf. is: %g\n", - begintimefit, integral); - - sprintf(buf, "test%d.xvg", nfitpnts); - fp = xvgropen(buf, "C(t) + Fit to C(t)", "Time (ps)", "C(t)", oenv); - fprintf(fp, "# parm0 = %g, parm1 = %g, parm2 = %g\n", - parm[0], parm[1], parm[2]); - for (j = 0; (j < nfitpnts); j++) - { - ttt = x0 ? x0[j] : dt*j; - fprintf(fp, "%10.5e %10.5e %10.5e\n", - ttt, c1[j], - lmcurves[eFitFn](ttt, parm)); - } - xvgrclose(fp); - } - } - if (fitparms) - { - for (i = 0; (i < nparm); i++) - { - fitparms[i] = parm[i]; - } - } - sfree(parm); - sfree(dparm); - } - - sfree(x); - sfree(y); - sfree(dy); - - return integral; -} - -/*! See description in header file. */ -real fit_acf(int ncorr, int fitfn, const output_env_t oenv, gmx_bool bVerbose, - real tbeginfit, real tendfit, real dt, real c1[], real *fit) -{ - double fitparm[3]; - double tStart, tail_corr, sum, sumtot = 0, c_start, ct_estimate; - real *sig; - int i, j, jmax, nf_int; - gmx_bool bPrint; - - bPrint = bVerbose || bDebugMode(); - - if (bPrint) - { - printf("COR:\n"); - } - - if (tendfit <= 0) - { - tendfit = ncorr*dt; - } - nf_int = min(ncorr, (int)(tendfit/dt)); - sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1); - - /* Estimate the correlation time for better fitting */ - ct_estimate = 0.5*c1[0]; - for (i = 1; (i < ncorr) && (c1[i] > 0); i++) - { - ct_estimate += c1[i]; - } - ct_estimate *= dt/c1[0]; - - if (bPrint) - { - printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n", - 0.0, dt*nf_int, sum); - printf("COR: Relaxation times are computed as fit to an exponential:\n"); - printf("COR: %s\n", longs_ffn[fitfn]); - printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n", tbeginfit, min(ncorr*dt, tendfit)); - } - - tStart = 0; - if (bPrint) - { - printf("COR:%11s%11s%11s%11s%11s%11s%11s\n", - "Fit from", "Integral", "Tail Value", "Sum (ps)", " a1 (ps)", - (nfp_ffn[fitfn] >= 2) ? " a2 ()" : "", - (nfp_ffn[fitfn] >= 3) ? " a3 (ps)" : ""); - } - - snew(sig, ncorr); - - if (tbeginfit > 0) - { - jmax = 3; - } - else - { - jmax = 1; - } - for (j = 0; ((j < jmax) && (tStart < tendfit) && (tStart < ncorr*dt)); j++) - { - /* Estimate the correlation time for better fitting */ - c_start = -1; - ct_estimate = 0; - for (i = 0; (i < ncorr) && (dt*i < tStart || c1[i] > 0); i++) - { - if (c_start < 0) - { - if (dt*i >= tStart) - { - c_start = c1[i]; - ct_estimate = 0.5*c1[i]; - } - } - else - { - ct_estimate += c1[i]; - } - } - if (c_start > 0) - { - ct_estimate *= dt/c_start; - } - else - { - /* The data is strange, so we need to choose somehting */ - ct_estimate = tendfit; - } - if (debug) - { - fprintf(debug, "tStart %g ct_estimate: %g\n", tStart, ct_estimate); - } - - if (fitfn == effnEXP3) - { - fitparm[0] = 0.002*ncorr*dt; - fitparm[1] = 0.95; - fitparm[2] = 0.2*ncorr*dt; - } - else - { - /* Good initial guess, this increases the probability of convergence */ - fitparm[0] = ct_estimate; - fitparm[1] = 1.0; - fitparm[2] = 1.0; - } - - /* Generate more or less appropriate sigma's */ - for (i = 0; i < ncorr; i++) - { - sig[i] = sqrt(ct_estimate+dt*i); - } - - nf_int = min(ncorr, (int)((tStart+1e-4)/dt)); - sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1); - tail_corr = do_lmfit(ncorr, c1, sig, dt, NULL, tStart, tendfit, oenv, - bDebugMode(), fitfn, fitparm, 0); - sumtot = sum+tail_corr; - if (fit && ((jmax == 1) || (j == 1))) - { - double mfp[3]; - for (i = 0; (i < 3); i++) - { - mfp[i] = fitparm[i]; - } - for (i = 0; (i < ncorr); i++) - { - fit[i] = lmcurves[fitfn](i*dt, mfp); - } - } - if (bPrint) - { - printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart, sum, tail_corr, sumtot); - for (i = 0; (i < nfp_ffn[fitfn]); i++) - { - printf(" %11.4e", fitparm[i]); - } - printf("\n"); - } - tStart += tbeginfit; - } - sfree(sig); - - return sumtot; -} diff --git a/src/gromacs/correlationfunctions/expfit.cpp b/src/gromacs/correlationfunctions/expfit.cpp new file mode 100644 index 0000000000..3b201bf2c9 --- /dev/null +++ b/src/gromacs/correlationfunctions/expfit.cpp @@ -0,0 +1,986 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2004, The GROMACS development team. + * Copyright (c) 2013,2014,2015, 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. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +/*! \internal + * \file + * \brief + * Implements routine for fitting a data set to a curve + * + * \author David van der Spoel + * \ingroup module_correlationfunctions + */ +#include "gmxpre.h" + +#include "expfit.h" + +#include +#include + +#include + +#include "external/lmfit/lmcurve.h" + +#include "gromacs/correlationfunctions/integrate.h" +#include "gromacs/fileio/xvgr.h" +#include "gromacs/legacyheaders/macros.h" +#include "gromacs/math/vec.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/futil.h" +#include "gromacs/utility/gmxassert.h" +#include "gromacs/utility/real.h" +#include "gromacs/utility/smalloc.h" + +/*! \brief Number of parameters for each fitting function */ +static const int nfp_ffn[effnNR] = { 0, 1, 2, 3, 5, 7, 9, 2, 4, 3, 6 }; + +/* +2 becuase parse_common_args wants leading and concluding NULL. + * We only allow exponential functions as choices on the command line, + * hence there are many more NULL field (which have to be at the end of + * the array). + */ +const char *s_ffn[effnNR+2] = { + NULL, "none", "exp", "aexp", "exp_exp", + "exp5", "exp7", "exp9", + NULL, NULL, NULL, NULL, NULL +}; + +/*! \brief Long description for each fitting function type */ +static const char *longs_ffn[effnNR] = { + "no fit", + "y = exp(-x/|a0|)", + "y = a1 exp(-x/|a0|)", + "y = a1 exp(-x/|a0|) + (1-a1) exp(-x/(|a2|)), a2 > a0 > 0", + "y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4, a3 >= a1", + "y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6, a5 >= a3 >= a1", + "y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6 exp(-x/|a7|) + a8, a7 >= a5 >= a3 >= a1", + "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a0), w = sqrt(1 - a1)", + "y = 1/2*(a0+a1) - 1/2*(a0-a1)*erf( (x-a2) /a3^2)", + "y = a1 *2*a0*((exp(-x/a0)-1)*(a0/x)+1)+(1-a1)*2*a2*((exp(-x/a2)-1)*(a2/x)+1)", + "y = (1-a0)*cos(x*a1)*exp(-(x/a2)^a3) + a0*exp(-(x/a4)^a5)" +}; + +int effnNparams(int effn) +{ + if ((0 <= effn) && (effn < effnNR)) + { + return nfp_ffn[effn]; + } + else + { + return -1; + } +} + +const char *effnDescription(int effn) +{ + if ((0 <= effn) && (effn < effnNR)) + { + return longs_ffn[effn]; + } + else + { + return NULL; + } +} + +int sffn2effn(const char **sffn) +{ + int eFitFn, i; + + eFitFn = 0; + for (i = 0; i < effnNR; i++) + { + if (sffn[i+1] && strcmp(sffn[0], sffn[i+1]) == 0) + { + eFitFn = i; + } + } + + return eFitFn; +} + +/*! \brief Compute exponential function A exp(-x/tau) */ +static double myexp(double x, double A, double tau) +{ + if ((A == 0) || (tau == 0)) + { + return 0; + } + return A*exp(-x/tau); +} + +/*! \brief Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2) */ +static double lmc_erffit (double x, const double *a) +{ + double erfarg; + double myerf; + + if (a[3] != 0) + { + erfarg = (x-a[2])/(a[3]*a[3]); + myerf = gmx_erfd(erfarg); + } + else + { + /* If a[3] == 0, a[3]^2 = 0 and the erfarg becomes +/- infinity */ + if (x < a[2]) + { + myerf = -1; + } + else + { + myerf = 1; + } + } + return 0.5*((a[0]+a[1]) - (a[0]-a[1])*myerf); +} + +/*! \brief Exponent function that prevents overflow */ +static double safe_exp(double x) +{ + double exp_max = 200; + double exp_min = -exp_max; + if (x <= exp_min) + { + return exp(exp_min); + } + else if (x >= exp_max) + { + return exp(exp_max); + } + else + { + return exp(x); + } +} + +/*! \brief Exponent minus 1 function that prevents overflow */ +static double safe_expm1(double x) +{ + double exp_max = 200; + double exp_min = -exp_max; + if (x <= exp_min) + { + return -1; + } + else if (x >= exp_max) + { + return exp(exp_max); + } + else + { + return gmx_expm1(x); + } +} + +/*! \brief Compute y = exp(-x/|a0|) */ +static double lmc_exp_one_parm(double x, const double *a) +{ + return safe_exp(-x/fabs(a[0])); +} + +/*! \brief Compute y = a1 exp(-x/|a0|) */ +static double lmc_exp_two_parm(double x, const double *a) +{ + return a[1]*safe_exp(-x/fabs(a[0])); +} + +/*! \brief Compute y = a1 exp(-x/|a0|) + (1-a1) exp(-x/|a2|) */ +static double lmc_exp_exp(double x, const double *a) +{ + double e1, e2; + + e1 = safe_exp(-x/fabs(a[0])); + e2 = safe_exp(-x/(fabs(a[0])+fabs(a[2]))); + return a[1]*e1 + (1-a[1])*e2; +} + +/*! \brief Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + a4 */ +static double lmc_exp_5_parm(double x, const double *a) +{ + double e1, e2; + + e1 = safe_exp(-x/fabs(a[1])); + e2 = safe_exp(-x/(fabs(a[1])+fabs(a[3]))); + return a[0]*e1 + a[2]*e2 + a[4]; +} + +/*! \brief Compute 7 parameter exponential function value. + * + * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + + * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6 + */ +static double lmc_exp_7_parm(double x, const double *a) +{ + double e1, e2, e3; + double fa1, fa3, fa5; + + fa1 = fabs(a[1]); + fa3 = fa1 + fabs(a[3]); + fa5 = fa3 + fabs(a[5]); + e1 = safe_exp(-x/fa1); + e2 = safe_exp(-x/fa3); + e3 = safe_exp(-x/fa5); + return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6]; +} + +/*! \brief Compute 9 parameter exponential function value. + * + * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + + * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6 exp(-x/(|a1|+|a3|+|a5|+|a7|)) + a8 + */ +static double lmc_exp_9_parm(double x, const double *a) +{ + double e1, e2, e3, e4; + double fa1, fa3, fa5, fa7; + + fa1 = fabs(a[1]); + fa3 = fa1 + fabs(a[3]); + fa5 = fa3 + fabs(a[5]); + fa7 = fa5 + fabs(a[7]); + + e1 = safe_exp(-x/fa1); + e2 = safe_exp(-x/fa3); + e3 = safe_exp(-x/fa5); + e4 = safe_exp(-x/fa7); + return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6]*e4 + a[8]; +} + +/*! \brief Compute y = (1-a0)*exp(-(x/|a2|)^|a3|)*cos(x*|a1|) + a0*exp(-(x/|a4|)^|a5|) */ +static double lmc_pres_6_parm(double x, const double *a) +{ + double term1, term2, term3; + double pow_max = 10; + + term3 = 0; + if ((a[4] != 0) && (a[0] != 0)) + { + double power = std::min(fabs(a[5]), pow_max); + term3 = a[0] * safe_exp(-pow((x/fabs(a[4])), power)); + } + + term1 = 1-a[0]; + term2 = 0; + if ((term1 != 0) && (a[2] != 0)) + { + double power = std::min(fabs(a[3]), pow_max); + term2 = safe_exp(-pow((x/fabs(a[2])), power)) * cos(x*fabs(a[1])); + } + + return term1*term2 + term3; +} + +/*! \brief Compute vac function */ +static double lmc_vac_2_parm(double x, const double *a) +{ + /* Fit to function + * + * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v) + * + * = exp(-v) (cosh(wv) + 1/w sinh(wv)) + * + * v = x/(2 a0) + * w = sqrt(1 - a1) + * + * For tranverse current autocorrelation functions: + * a0 = tau + * a1 = 4 tau (eta/rho) k^2 + * + */ + + double y, v, det, omega, wv, em, ec, es; + double wv_max = 100; + + v = x/(2*fabs(a[0])); + det = 1 - a[1]; + em = safe_exp(-v); + if (det != 0) + { + omega = sqrt(fabs(det)); + wv = std::min(omega*v, wv_max); + + if (det > 0) + { + ec = em*0.5*(safe_exp(wv)+safe_exp(-wv)); + es = em*0.5*(safe_exp(wv)-safe_exp(-wv))/omega; + } + else + { + ec = em*cos(wv); + es = em*sin(wv)/omega; + } + y = ec + es; + } + else + { + y = (1+v)*em; + } + return y; +} + +/*! \brief Compute error estimate */ +static double lmc_errest_3_parm(double x, const double *a) +{ + double e1, e2, v1, v2; + double fa0 = fabs(a[0]); + double fa1; + double fa2 = fa0+fabs(a[2]); + + if (a[0] != 0) + { + e1 = safe_expm1(-x/fa0); + } + else + { + e1 = 0; + } + if (a[2] != 0) + { + e2 = safe_expm1(-x/fa2); + } + else + { + e2 = 0; + } + + if (x > 0) + { + v1 = 2*fa0*(e1*fa0/x + 1); + v2 = 2*fa2*(e2*fa2/x + 1); + /* We need 0 <= a1 <= 1 */ + fa1 = std::min(1.0, std::max(0.0, a[1])); + + return fa1*v1 + (1-fa1)*v2; + } + else + { + return 0; + } +} + +/*! \brief function type for passing to fitting routine */ +typedef double (*t_lmcurve)(double x, const double *a); + +/*! \brief array of fitting functions corresponding to the pre-defined types */ +t_lmcurve lmcurves[effnNR+1] = { + lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm, + lmc_exp_exp, lmc_exp_5_parm, lmc_exp_7_parm, + lmc_exp_9_parm, + lmc_vac_2_parm, lmc_erffit, lmc_errest_3_parm, lmc_pres_6_parm +}; + +double fit_function(const int eFitFn, const double parm[], const double x) +{ + if ((eFitFn < 0) || (eFitFn >= effnNR)) + { + fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n", + eFitFn, effnNR-1); + return 0.0; + } + return lmcurves[eFitFn](x, parm); +} + +/*! \brief lmfit_exp supports fitting of different functions + * + * This routine calls the Levenberg-Marquardt non-linear fitting + * routine for fitting a data set with errors to a target function. + * Fitting routines included in gromacs in src/external/lmfit. + */ +static gmx_bool lmfit_exp(int nfit, + const double x[], + const double y[], + const double dy[], + double parm[], + gmx_bool bVerbose, + int eFitFn, + int nfix) +{ + double chisq, ochisq; + gmx_bool bCont; + int j; + int maxiter = 100; + lm_control_struct control; + lm_status_struct *status; + int nparam = effnNparams(eFitFn); + int p2; + gmx_bool bSkipLast; + + if ((eFitFn < 0) || (eFitFn >= effnNR)) + { + fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n", + eFitFn, effnNR-1); + return FALSE; + } + /* Using default control structure for double precision fitting that + * comes with the lmfit package (i.e. from the include file). + */ + control = lm_control_double; + control.verbosity = (bVerbose ? 1 : 0); + control.n_maxpri = 0; + control.m_maxpri = 0; + + snew(status, 1); + /* Initial params */ + chisq = 1e12; + j = 0; + if (bVerbose) + { + printf("%4s %10s Parameters\n", "Step", "chi^2"); + } + /* Check whether we have to skip some params */ + if (nfix > 0) + { + do + { + p2 = 1 << (nparam-1); + bSkipLast = ((p2 & nfix) == p2); + if (bSkipLast) + { + nparam--; + nfix -= p2; + } + } + while ((nparam > 0) && (bSkipLast)); + if (bVerbose) + { + printf("Using %d out of %d parameters\n", nparam, effnNparams(eFitFn)); + } + } + do + { + ochisq = chisq; + lmcurve(nparam, parm, nfit, x, y, dy, + lmcurves[eFitFn], &control, status); + chisq = sqr(status->fnorm); + if (bVerbose) + { + printf("status: fnorm = %g, nfev = %d, userbreak = %d\noutcome = %s\n", + status->fnorm, status->nfev, status->userbreak, + lm_infmsg[status->outcome]); + } + if (bVerbose) + { + int mmm; + printf("%4d %8g", j, chisq); + for (mmm = 0; (mmm < effnNparams(eFitFn)); mmm++) + { + printf(" %8g", parm[mmm]); + } + printf("\n"); + } + j++; + bCont = (fabs(ochisq - chisq) > fabs(control.ftol*chisq)); + } + while (bCont && (j < maxiter)); + + sfree(status); + + return TRUE; +} + +/*! \brief Ensure the fitting parameters are well-behaved. + * + * In order to make sure that fitting behaves according to the + * constraint that time constants are positive and increasing + * we enforce this here by setting all time constants to their + * absolute value and by adding e.g. |a_0| to |a_2|. This is + * done in conjunction with the fitting functions themselves. + * When there are multiple time constants we make sure that + * the successive time constants are at least double the + * previous ones and during fitting we enforce the they remain larger. + * It may very well help the convergence of the fitting routine. + */ +static void initiate_fit_params(int eFitFn, + double params[]) +{ + int i, nparm; + + nparm = effnNparams(eFitFn); + + switch (eFitFn) + { + case effnVAC: + GMX_ASSERT(params[0] >= 0, "parameters should be >= 0"); + break; + case effnEXP1: + case effnEXP2: + case effnEXPEXP: + GMX_ASSERT(params[0] >= 0, "parameters should be >= 0"); + if (nparm > 2) + { + GMX_ASSERT(params[2] >= 0, "parameters should be >= 0"); + /* In order to maintain params[2] >= params[0] in the final + * result, we fit the difference between the two, that is + * params[2]-params[0] and in the add add in params[0] + * again. + */ + params[2] = std::max(fabs(params[2])-params[0], params[0]); + } + break; + case effnEXP5: + case effnEXP7: + case effnEXP9: + GMX_ASSERT(params[1] >= 0, "parameters should be >= 0"); + params[1] = fabs(params[1]); + if (nparm > 3) + { + GMX_ASSERT(params[3] >= 0, "parameters should be >= 0"); + /* See comment under effnEXPEXP */ + params[3] = std::max(fabs(params[3])-params[1], params[1]); + if (nparm > 5) + { + GMX_ASSERT(params[5] >= 0, "parameters should be >= 0"); + /* See comment under effnEXPEXP */ + params[5] = std::max(fabs(params[5])-params[3], params[3]); + if (nparm > 7) + { + GMX_ASSERT(params[7] >= 0, "parameters should be >= 0"); + /* See comment under effnEXPEXP */ + params[7] = std::max(fabs(params[7])-params[5], params[5]); + } + } + } + break; + case effnERREST: + GMX_ASSERT(params[0] >= 0, "parameters should be >= 0"); + GMX_ASSERT(params[1] >= 0 && params[1] <= 1, "parameter 1 should in 0 .. 1"); + GMX_ASSERT(params[2] >= 0, "parameters should be >= 0"); + /* See comment under effnEXPEXP */ + params[2] = fabs(params[2])-params[0]; + break; + case effnPRES: + for (i = 1; (i < nparm); i++) + { + GMX_ASSERT(params[i] >= 0, "parameters should be >= 0"); + } + break; + default: + break; + } +} + +/*! \brief Process the fitting parameters to get output parameters. + * + * See comment at the previous function. + */ +static void extract_fit_params(int eFitFn, + double params[]) +{ + int i, nparm; + + nparm = effnNparams(eFitFn); + + switch (eFitFn) + { + case effnVAC: + params[0] = fabs(params[0]); + break; + case effnEXP1: + case effnEXP2: + case effnEXPEXP: + params[0] = fabs(params[0]); + if (nparm > 2) + { + /* Back conversion of parameters from the fitted difference + * to the absolute value. + */ + params[2] = fabs(params[2])+params[0]; + } + break; + case effnEXP5: + case effnEXP7: + case effnEXP9: + params[1] = fabs(params[1]); + if (nparm > 3) + { + /* See comment under effnEXPEXP */ + params[3] = fabs(params[3])+params[1]; + if (nparm > 5) + { + /* See comment under effnEXPEXP */ + params[5] = fabs(params[5])+params[3]; + if (nparm > 7) + { + /* See comment under effnEXPEXP */ + params[7] = fabs(params[7])+params[5]; + } + } + } + break; + case effnERREST: + params[0] = fabs(params[0]); + if (params[1] < 0) + { + params[1] = 0; + } + else if (params[1] > 1) + { + params[1] = 1; + } + /* See comment under effnEXPEXP */ + params[2] = params[0]+fabs(params[2]); + break; + case effnPRES: + for (i = 1; (i < nparm); i++) + { + params[i] = fabs(params[i]); + } + break; + default: + break; + } +} + +/*! \brief Print chi-squared value and the parameters */ +static void print_chi2_params(FILE *fp, + const int eFitFn, + const double fitparms[], + const char *label, + const int nfitpnts, + const double x[], + const double y[]) +{ + int i; + double chi2 = 0; + + for (i = 0; (i < nfitpnts); i++) + { + double yfit = lmcurves[eFitFn](x[i], fitparms); + chi2 += sqr(y[i] - yfit); + } + fprintf(fp, "There are %d data points, %d parameters, %s chi2 = %g\nparams:", + nfitpnts, effnNparams(eFitFn), label, chi2); + for (i = 0; (i < effnNparams(eFitFn)); i++) + { + fprintf(fp, " %10g", fitparms[i]); + } + fprintf(fp, "\n"); +} + +/*! \brief See description in header file. */ +real do_lmfit(int ndata, real c1[], real sig[], real dt, real x0[], + real begintimefit, real endtimefit, const output_env_t oenv, + gmx_bool bVerbose, int eFitFn, double fitparms[], int fix, + const char *fn_fitted) +{ + FILE *fp; + int i, j, nfitpnts; + double integral, ttt; + double *x, *y, *dy; + + if (0 != fix) + { + fprintf(stderr, "Using fixed parameters in curve fitting is temporarily not working.\n"); + } + if (debug) + { + fprintf(debug, "There are %d points to fit %d vars!\n", ndata, effnNparams(eFitFn)); + fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n", + eFitFn, begintimefit, endtimefit, dt); + } + + snew(x, ndata); + snew(y, ndata); + snew(dy, ndata); + j = 0; + for (i = 0; i < ndata; i++) + { + ttt = x0 ? x0[i] : dt*i; + if (ttt >= begintimefit && ttt <= endtimefit) + { + x[j] = ttt; + y[j] = c1[i]; + if (NULL == sig) + { + // No weighting if all values are divided by 1. + dy[j] = 1; + } + else + { + dy[j] = std::max(1.0e-7, (double)sig[i]); + } + if (debug) + { + fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy=%g, ttt=%g\n", + j, i, x[j], y[j], dy[j], ttt); + } + j++; + } + } + nfitpnts = j; + integral = 0; + if (nfitpnts < effnNparams(eFitFn)) + { + fprintf(stderr, "Not enough (%d) data points for fitting, dt = %g!\n", + nfitpnts, dt); + } + else + { + gmx_bool bSuccess; + + if (bVerbose) + { + print_chi2_params(stdout, eFitFn, fitparms, "initial", nfitpnts, x, y); + } + initiate_fit_params(eFitFn, fitparms); + + bSuccess = lmfit_exp(nfitpnts, x, y, dy, fitparms, bVerbose, eFitFn, fix); + extract_fit_params(eFitFn, fitparms); + + if (!bSuccess) + { + fprintf(stderr, "Fit failed!\n"); + } + else + { + if (bVerbose) + { + print_chi2_params(stdout, eFitFn, fitparms, "final", nfitpnts, x, y); + } + switch (eFitFn) + { + case effnEXP1: + integral = fitparms[0]*myexp(begintimefit, 1, fitparms[0]); + break; + case effnEXP2: + integral = fitparms[0]*myexp(begintimefit, fitparms[1], fitparms[0]); + break; + case effnEXPEXP: + integral = (fitparms[0]*myexp(begintimefit, fitparms[1], fitparms[0]) + + fitparms[2]*myexp(begintimefit, 1-fitparms[1], fitparms[2])); + break; + case effnEXP5: + case effnEXP7: + case effnEXP9: + integral = 0; + for (i = 0; (i < (effnNparams(eFitFn)-1)/2); i++) + { + integral += fitparms[2*i]*myexp(begintimefit, fitparms[2*i+1], fitparms[2*i]); + } + break; + default: + /* Do numerical integration */ + integral = 0; + for (i = 0; (i < nfitpnts-1); i++) + { + double y0 = lmcurves[eFitFn](x[i], fitparms); + double y1 = lmcurves[eFitFn](x[i+1], fitparms); + integral += (x[i+1]-x[i])*(y1+y0)*0.5; + } + break; + } + + if (bVerbose) + { + printf("FIT: Integral of fitted function: %g\n", integral); + if ((effnEXP5 == eFitFn) || (effnEXP7 == eFitFn) || (effnEXP9 == eFitFn)) + { + printf("FIT: Note that the constant term is not taken into account when computing integral.\n"); + } + } + /* Generate debug output */ + if (NULL != fn_fitted) + { + fp = xvgropen(fn_fitted, "Data + Fit", "Time (ps)", + "Data (t)", oenv); + for (i = 0; (i < effnNparams(eFitFn)); i++) + { + fprintf(fp, "# fitparms[%d] = %g\n", i, fitparms[i]); + } + for (j = 0; (j < nfitpnts); j++) + { + real ttt = x0 ? x0[i] : dt*j; + fprintf(fp, "%10.5e %10.5e %10.5e\n", + x[j], y[j], lmcurves[eFitFn](ttt, fitparms)); + } + xvgrclose(fp); + } + } + } + + sfree(x); + sfree(y); + sfree(dy); + + return integral; +} + +/*! See description in header file. */ +real fit_acf(int ncorr, int fitfn, const output_env_t oenv, gmx_bool bVerbose, + real tbeginfit, real tendfit, real dt, real c1[], real *fit) +{ + double fitparm[3]; + double tStart, tail_corr, sum, sumtot = 0, c_start, ct_estimate; + real *sig; + int i, j, jmax, nf_int; + gmx_bool bPrint; + + bPrint = bVerbose || bDebugMode(); + + if (bPrint) + { + printf("COR:\n"); + } + + if (tendfit <= 0) + { + tendfit = ncorr*dt; + } + nf_int = std::min(ncorr, (int)(tendfit/dt)); + sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1); + + /* Estimate the correlation time for better fitting */ + ct_estimate = 0.5*c1[0]; + for (i = 1; (i < ncorr) && (c1[i] > 0); i++) + { + ct_estimate += c1[i]; + } + ct_estimate *= dt/c1[0]; + + if (bPrint) + { + printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n", + 0.0, dt*nf_int, sum); + printf("COR: Relaxation times are computed as fit to an exponential:\n"); + printf("COR: %s\n", effnDescription(fitfn)); + printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n", tbeginfit, std::min(ncorr*dt, tendfit)); + } + + tStart = 0; + if (bPrint) + { + printf("COR:%11s%11s%11s%11s%11s%11s%11s\n", + "Fit from", "Integral", "Tail Value", "Sum (ps)", " a1 (ps)", + (effnNparams(fitfn) >= 2) ? " a2 ()" : "", + (effnNparams(fitfn) >= 3) ? " a3 (ps)" : ""); + } + + snew(sig, ncorr); + + if (tbeginfit > 0) + { + jmax = 3; + } + else + { + jmax = 1; + } + for (j = 0; ((j < jmax) && (tStart < tendfit) && (tStart < ncorr*dt)); j++) + { + /* Estimate the correlation time for better fitting */ + c_start = -1; + ct_estimate = 0; + for (i = 0; (i < ncorr) && (dt*i < tStart || c1[i] > 0); i++) + { + if (c_start < 0) + { + if (dt*i >= tStart) + { + c_start = c1[i]; + ct_estimate = 0.5*c1[i]; + } + } + else + { + ct_estimate += c1[i]; + } + } + if (c_start > 0) + { + ct_estimate *= dt/c_start; + } + else + { + /* The data is strange, so we need to choose somehting */ + ct_estimate = tendfit; + } + if (debug) + { + fprintf(debug, "tStart %g ct_estimate: %g\n", tStart, ct_estimate); + } + + if (fitfn == effnEXPEXP) + { + fitparm[0] = 0.002*ncorr*dt; + fitparm[1] = 0.95; + fitparm[2] = 0.2*ncorr*dt; + } + else + { + /* Good initial guess, this increases the probability of convergence */ + fitparm[0] = ct_estimate; + fitparm[1] = 1.0; + fitparm[2] = 1.0; + } + + /* Generate more or less appropriate sigma's */ + for (i = 0; i < ncorr; i++) + { + sig[i] = sqrt(ct_estimate+dt*i); + } + + nf_int = std::min(ncorr, (int)((tStart+1e-4)/dt)); + sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1); + tail_corr = do_lmfit(ncorr, c1, sig, dt, NULL, tStart, tendfit, oenv, + bDebugMode(), fitfn, fitparm, 0, NULL); + sumtot = sum+tail_corr; + if (fit && ((jmax == 1) || (j == 1))) + { + double mfp[3]; + for (i = 0; (i < 3); i++) + { + mfp[i] = fitparm[i]; + } + for (i = 0; (i < ncorr); i++) + { + fit[i] = lmcurves[fitfn](i*dt, mfp); + } + } + if (bPrint) + { + printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart, sum, tail_corr, sumtot); + for (i = 0; (i < effnNparams(fitfn)); i++) + { + printf(" %11.4e", fitparm[i]); + } + printf("\n"); + } + tStart += tbeginfit; + } + sfree(sig); + + return sumtot; +} diff --git a/src/gromacs/correlationfunctions/expfit.h b/src/gromacs/correlationfunctions/expfit.h index ee75f8d711..e9cb14938d 100644 --- a/src/gromacs/correlationfunctions/expfit.h +++ b/src/gromacs/correlationfunctions/expfit.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2014, by the GROMACS development team, led by + * Copyright (c) 2014,2015, 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. @@ -55,8 +55,9 @@ extern "C" { * Enum to select fitting functions */ enum { - effnNONE, effnEXP1, effnEXP2, effnEXP3, effnVAC, - effnEXP5, effnEXP7, effnEXP9, effnERF, effnERREST, effnPRES, effnNR + effnNONE, effnEXP1, effnEXP2, effnEXPEXP, + effnEXP5, effnEXP7, effnEXP9, + effnVAC, effnERF, effnERREST, effnPRES, effnNR }; /*! \brief @@ -94,7 +95,7 @@ int sffn2effn(const char **sffn); * \param[in] x The value of x * \return the value of the fit */ -double fit_function(int eFitFn, double *parm, double x); +double fit_function(const int eFitFn, const double parm[], const double x); /*! \brief * Use Levenberg-Marquardt method to fit to a nfitparm parameter exponential @@ -103,7 +104,7 @@ double fit_function(int eFitFn, double *parm, double x); * If x == NULL, the timestep dt will be used to create a time axis. * \param[in] ndata Number of data points * \param[in] c1 The data points - * \param[in] sig The standard deviation in the points + * \param[in] sig The standard deviation in the points (can be NULL) * \param[in] dt The time step * \param[in] x The X-axis (may be NULL, see above) * \param[in] begintimefit Starting time for fitting @@ -111,14 +112,22 @@ double fit_function(int eFitFn, double *parm, double x); * \param[in] oenv Output formatting information * \param[in] bVerbose Should the routine write to console? * \param[in] eFitFn Fitting function (0 .. effnNR) - * \param[out] fitparms[] + * \param[inout] fitparms[] Fitting parameters, see printed manual for a + * detailed description. Note that in all implemented cases the parameters + * corresponding to time constants will be generated with increasing values. + * Such input parameters should therefore be provided in increasing order. + * If this is not the case or if subsequent time parameters differ by less than + * a factor of 2, they will be modified to ensure tau_i+1 >= 2 tau_i. * \param[in] fix Constrains fit parameter i at it's starting value, when the i'th bit - * of fix is set. + * of fix is set. This works only when the N last parameters are fixed + * but not when a parameter somewhere in the middle needs to be fixed. + * \param[in] fn_fitted If not NULL file to print the data and fitted curve to * \return integral. */ real do_lmfit(int ndata, real c1[], real sig[], real dt, real *x, real begintimefit, real endtimefit, const output_env_t oenv, - gmx_bool bVerbose, int eFitFn, double fitparms[], int fix); + gmx_bool bVerbose, int eFitFn, double fitparms[], int fix, + const char *fn_fitted); /*! \brief * Fit an autocorrelation function to a pre-defined functional form @@ -132,7 +141,7 @@ real do_lmfit(int ndata, real c1[], real sig[], real dt, real *x, * \param[in] tendfit Ending time for fitting * \param[in] dt The time step * \param[in] c1 The data points - * \param[out] fit The fitting parameters + * \param[inout] fit The fitting parameters * \return the integral over the autocorrelation function? */ real fit_acf(int ncorr, int fitfn, const output_env_t oenv, gmx_bool bVerbose, diff --git a/src/gromacs/correlationfunctions/tests/expfit.cpp b/src/gromacs/correlationfunctions/tests/expfit.cpp old mode 100755 new mode 100644 index c9f518b392..74fb2ebd8e --- a/src/gromacs/correlationfunctions/tests/expfit.cpp +++ b/src/gromacs/correlationfunctions/tests/expfit.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2014, by the GROMACS development team, led by + * Copyright (c) 2014,2015, 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. @@ -49,34 +49,34 @@ #include #include "gromacs/fileio/xvgr.h" +#include "gromacs/legacyheaders/oenv.h" #include "gromacs/utility/smalloc.h" #include "testutils/refdata.h" #include "testutils/testasserts.h" #include "testutils/testfilemanager.h" -//! Number of data files for testing. -#define expTestNrTypes 3 - namespace gmx { namespace { +class ExpfitData +{ + public: + int nrLines_; + std::vector x_, y_; + real startTime_, endTime_, dt_; +}; + class ExpfitTest : public ::testing::Test { protected: - static int nrLines_; - static std::vector values_[expTestNrTypes]; - static int nrColumns_; - static std::vector standardDev_; - static real startTime_; - static real endTime_; - static real timeDeriv_; - test::TestReferenceData refData_; - test::TestReferenceChecker checker_; + static std::vector data_; + test::TestReferenceData refData_; + test::TestReferenceChecker checker_; ExpfitTest( ) : checker_(refData_.rootChecker()) { @@ -85,40 +85,34 @@ class ExpfitTest : public ::testing::Test // Static initiation, only run once every test. static void SetUpTestCase() { - double ** tempValues; - std::string fileName[expTestNrTypes]; - fileName[0] = test::TestFileManager::getInputFilePath("testINVEXP.xvg"); - fileName[1] = test::TestFileManager::getInputFilePath("testPRES.xvg"); - fileName[2] = test::TestFileManager::getInputFilePath("testEXP.xvg"); - for (int i = 0; i < expTestNrTypes; i++) + double ** tempValues = NULL; + std::vector fileName; + fileName.push_back(test::TestFileManager::getInputFilePath("testINVEXP.xvg")); + fileName.push_back(test::TestFileManager::getInputFilePath("testPRES.xvg")); + fileName.push_back(test::TestFileManager::getInputFilePath("testINVEXP79.xvg")); + fileName.push_back(test::TestFileManager::getInputFilePath("testERF.xvg")); + fileName.push_back(test::TestFileManager::getInputFilePath("testERREST.xvg")); + for (std::vector::iterator i = fileName.begin(); i < fileName.end(); ++i) { - const char * name = fileName[i].c_str(); - // TODO: this assumes all files have the same length. - nrLines_ = read_xvg(name, &tempValues, &nrColumns_); - - // Generating standard deviation - if (i == 0) + const char * name = i->c_str(); + int nrColumns; + ExpfitData ed; + ed.nrLines_ = read_xvg(name, &tempValues, &nrColumns); + ed.dt_ = tempValues[0][1] - tempValues[0][0]; + ed.startTime_ = tempValues[0][0]; + ed.endTime_ = tempValues[0][ed.nrLines_-1]; + for (int j = 0; j < ed.nrLines_; j++) { - double fac = 1.0/nrLines_; - for (int j = 0; j < nrLines_; j++) - { - standardDev_.push_back(fac); - } - timeDeriv_ = tempValues[0][1] - tempValues[0][0]; - startTime_ = tempValues[0][0]; - endTime_ = tempValues[0][nrLines_-1]; - } - - for (int j = 0; j < nrLines_; j++) - { - values_[i].push_back((real)tempValues[1][j]); + ed.x_.push_back((real)tempValues[0][j]); + ed.y_.push_back((real)tempValues[1][j]); } + data_.push_back(ed); // Free memory that was allocated in read_xvg - for (int i = 0; i < nrColumns_; i++) + for (int j = 0; j < nrColumns; j++) { - sfree(tempValues[i]); - tempValues[i] = NULL; + sfree(tempValues[j]); + tempValues[j] = NULL; } sfree(tempValues); tempValues = NULL; @@ -129,13 +123,26 @@ class ExpfitTest : public ::testing::Test { } - void test(int type, double result[], double tolerance, int testType) + void test(int type, double result[], double tolerance, + unsigned int testType) { - int nfitparm = effnNparams(type); - - do_lmfit(nrLines_, &values_[testType][0], &standardDev_[0], timeDeriv_, - NULL, startTime_, endTime_, NULL, false, type, result, 0); + int nfitparm = effnNparams(type); + output_env_t oenv; + if (testType >= data_.size()) + { + GMX_THROW(InvalidInputError("testType out of range")); + } + output_env_init_default(&oenv); + do_lmfit(data_[testType].nrLines_, + &(data_[testType].y_[0]), + NULL, + data_[testType].dt_, + &(data_[testType].x_[0]), + data_[testType].startTime_, + data_[testType].endTime_, + oenv, false, type, result, 0, NULL); + output_env_done(oenv); checker_.setDefaultTolerance(test::relativeToleranceAsFloatingPoint(1, tolerance)); checker_.checkSequenceArray(nfitparm, result, "result"); } @@ -143,63 +150,55 @@ class ExpfitTest : public ::testing::Test //static var -int ExpfitTest::nrLines_; -//cppcheck-suppress arrayIndexOutOfBounds fixed in 1.68-dev -std::vector ExpfitTest::values_[expTestNrTypes]; -int ExpfitTest::nrColumns_; -std::vector ExpfitTest::standardDev_; -real ExpfitTest::startTime_; -real ExpfitTest::endTime_; -real ExpfitTest::timeDeriv_; +std::vector ExpfitTest::data_; TEST_F (ExpfitTest, EffnEXP1) { double param[] = {25}; - test(effnEXP1, param, 1e-6, 0); + test(effnEXP1, param, 1e-5, 0); } TEST_F (ExpfitTest, EffnEXP2) { double param[] = {35, 0.5}; - test(effnEXP2, param, 1e-6, 0); + test(effnEXP2, param, 3e-5, 0); } -TEST_F (ExpfitTest, EffnEXP3) { - double param[] = {45, 0.5, 5}; - test(effnEXP3, param, 1e-4, 0); +TEST_F (ExpfitTest, EffnEXPEXP) { + double param[] = {5, 0.5, 45}; + test(effnEXPEXP, param, 1e-2, 0); } TEST_F (ExpfitTest, EffnEXP5) { double param[] = {0.5, 5, 0.5, 50, 0.002}; - test(effnEXP5, param, 1e-4, 0); + test(effnEXP5, param, 1e-2, 2); } TEST_F (ExpfitTest, EffnEXP7) { - double param[] = {0.5, 5, -0.02, 0.5, 0.5, 50, -0.002}; - test(effnEXP7, param, 1e-4, 0); + double param[] = {0.1, 2, 0.5, 30, 0.3, 50, -0.002}; + test(effnEXP7, param, 1e-2, 2); } TEST_F (ExpfitTest, EffnEXP9) { - double param[] = {2, 1200, -1, 300, 0.7, 70, 0.5, 6, -0.5}; - test(effnEXP9, param, 4e-2, 0); + double param[] = {0.4, 5, 0.2, 30, 0.1, 70, 0.2, 200, -0.05}; + test(effnEXP9, param, 4e-2, 2); } TEST_F (ExpfitTest, EffnERF) { - double param[] = {0.5, 0.5, 0.5, 1}; - test(effnERF, param, 1e-2, 0); + double param[] = {80, 120, 180, 5}; + test(effnERF, param, 1e-1, 3); } TEST_F (ExpfitTest, EffnERREST) { - double param[] = {0.5, 0.7, 0.3}; - test(effnERREST, param, 1e-4, 2); + double param[] = {1, 0.9, 100}; + test(effnERREST, param, 5e-3, 4); } TEST_F (ExpfitTest, EffnVAC) { - double param[] = {0.5, 0.05}; - test(effnVAC, param, 1e-4, 0); + double param[] = {30, 0.0}; + test(effnVAC, param, 0.05, 0); } -TEST_F (ExpfitTest, DISABLED_EffnPRES) { - //TODO: This test is prodocues NaNs and INFs. Fix and then reactivate. - double param[] = {0, 10, 4, 1, 0.5, 1}; +TEST_F (ExpfitTest, EffnPRES) { + double param[] = {0.6, 10, 7, 1, 0.25, 2}; test(effnPRES, param, 1e-4, 1); } diff --git a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacCos.xml b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacCos.xml index 34cba4032b..41316e618e 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacCos.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacCos.xml @@ -14,7 +14,7 @@ 0.98145622 0.97775453 0.97388524 - 0.96959651 + 0.96959645 0.96521682 0.96071005 0.95656222 @@ -37,7 +37,7 @@ 0.92004192 0.92134565 0.92249507 - 0.92415452 + 0.92415446 0.92598462 0.92835307 0.93056071 @@ -48,8 +48,8 @@ 0.9445464 0.94751155 0.95099658 - 0.95437491 - 0.95766789 + 0.95437485 + 0.95766783 0.96085715 0.96405703 0.96722084 @@ -75,19 +75,19 @@ 0.98345906 0.98124063 0.97880316 - 0.97668463 + 0.97668457 0.97459376 0.97177458 0.9690451 0.96629435 - 0.96378684 + 0.96378678 0.96081358 0.95825708 0.95578271 0.95307195 - 0.95074266 + 0.9507426 0.94830912 - 0.94601887 + 0.94601882 0.94421709 0.94242907 0.94094354 @@ -106,9 +106,9 @@ 0.94294119 0.94477904 0.94679701 - 0.94860345 + 0.94860339 0.95079654 - 0.95264018 + 0.95264012 0.95462155 0.95684212 0.95912236 @@ -133,14 +133,14 @@ 0.97891223 0.97828108 0.97720605 - 0.97623944 + 0.97623938 0.97527778 0.97370166 0.97211659 0.97065586 0.96862853 0.96749282 - 0.96519589 + 0.96519583 0.96352464 0.96173507 0.9598251 @@ -155,7 +155,7 @@ 0.94696826 0.94605809 0.94519788 - 0.94469965 + 0.94469959 0.94453776 0.94388711 0.94351792 @@ -164,21 +164,21 @@ 0.94424397 0.94460499 0.94499773 - 0.94563341 + 0.94563335 0.94637001 0.94729853 0.94789016 - 0.9492957 + 0.94929564 0.95013034 0.95128489 - 0.95289272 + 0.95289266 0.9538911 0.95517474 0.95651346 0.95739591 0.95905453 - 0.95995909 - 0.96129471 + 0.95995903 + 0.96129465 0.96225822 0.96320486 0.96370006 @@ -194,7 +194,7 @@ 0.96682173 0.96666759 0.96608764 - 0.96583378 + 0.96583372 0.96554524 0.9645654 0.96389264 @@ -228,7 +228,7 @@ 0.94219875 0.94211328 0.94299912 - 0.94275576 + 0.9427557 0.94383633 0.94452924 0.94524854 @@ -505,5 +505,5 @@ 0.00094952446 0.00096428557 - 240.60057076136582 + 240.60056956927292 diff --git a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacNormal.xml b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacNormal.xml index 415cf1f556..b61d14bc34 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacNormal.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacNormal.xml @@ -12,26 +12,26 @@ 0.70085573 0.63575351 0.56621635 - 0.49021092 - 0.41346616 + 0.49021095 + 0.41346619 0.32695487 0.23947404 0.15114491 - 0.068911292 - -0.018343732 + 0.068911299 + -0.018343734 -0.10222372 -0.17851034 - -0.2556752 + -0.25567523 -0.33233142 -0.40906069 - -0.47129509 + -0.47129512 -0.52546626 -0.577564 -0.62862134 -0.66793829 -0.70031446 -0.71649045 - -0.73871797 + -0.73871803 -0.74570811 -0.75775838 -0.75111836 @@ -39,7 +39,7 @@ -0.71327382 -0.68468153 -0.65119308 - -0.60643828 + -0.60643834 -0.56428039 -0.52878588 -0.46647781 @@ -47,75 +47,75 @@ -0.34746397 -0.28559166 -0.22607948 - -0.15521984 + -0.15521985 -0.087527938 - -0.021420239 + -0.02142024 0.042880654 0.10628833 - 0.169659 + 0.16965902 0.23284768 - 0.29287171 + 0.29287174 0.33996987 0.3846795 0.4326345 0.47152123 0.50605392 - 0.53039247 + 0.53039253 0.552634 0.57193202 0.57857972 0.58289105 - 0.5867179 + 0.58671796 0.57858908 0.57083881 0.54604369 0.52723479 0.50264382 0.4660213 - 0.42628899 + 0.42628902 0.38307241 - 0.33706102 + 0.33706105 0.29605037 0.25652972 - 0.20322086 - 0.15179904 + 0.20322087 + 0.15179905 0.1005994 0.053715136 - -0.0024519849 + -0.0024519851 -0.051181242 - -0.098006107 + -0.098006114 -0.15021549 - -0.19438037 + -0.19438039 -0.2419212 - -0.28681141 + -0.28681144 -0.32269096 -0.35785615 - -0.38878986 + -0.38878989 -0.41733515 - -0.44138661 - -0.44811302 + -0.44138664 + -0.44811305 -0.46269119 -0.46645734 - -0.47144547 - -0.47026828 - -0.46132591 + -0.4714455 + -0.47026831 + -0.46132594 -0.44804862 -0.43905067 -0.41386977 - -0.38871932 - -0.36298344 + -0.38871935 + -0.36298347 -0.32728449 -0.28807318 -0.25138158 - -0.20809993 + -0.20809995 -0.17095725 -0.13121751 -0.086369112 -0.040557165 0.0057886061 - 0.050152354 - 0.082554683 - 0.12262377 + 0.050152358 + 0.082554691 + 0.12262378 0.15991174 0.19901413 0.22433081 @@ -130,82 +130,82 @@ 0.37430125 0.37234989 0.36795545 - 0.35466814 - 0.34249428 + 0.35466817 + 0.34249431 0.32283241 0.30485228 0.28708842 - 0.25817123 - 0.22833388 + 0.25817126 + 0.22833389 0.20199309 - 0.16451561 - 0.14491151 + 0.16451563 + 0.14491153 0.10192882 - 0.071294688 - 0.03886297 - 0.0038909942 + 0.071294695 + 0.038862973 + 0.0038909945 -0.020299936 -0.055376183 -0.092676423 -0.12172776 -0.15591128 - -0.17435807 + -0.17435808 -0.19855703 -0.21759604 - -0.2342567 + -0.23425671 -0.2501609 - -0.26639262 + -0.26639265 -0.27480033 -0.2769354 -0.28888217 -0.29434624 - -0.29110822 + -0.29110825 -0.28738263 -0.27644756 -0.26865044 - -0.26017156 - -0.24633232 + -0.26017159 + -0.24633233 -0.22993366 - -0.21036342 + -0.21036343 -0.19680807 -0.16794714 - -0.14982072 + -0.14982073 -0.1253354 - -0.091874786 + -0.091874793 -0.070614971 -0.043464452 -0.01496042 0.0054097059 - 0.039536409 + 0.039536413 0.058995754 0.087431222 - 0.10856776 - 0.12864889 + 0.10856777 + 0.12864891 0.14001396 0.15372792 0.16890197 0.1832682 - 0.19547278 + 0.19547279 0.20601526 - 0.22224233 + 0.22224234 0.22276327 - 0.22658545 + 0.22658546 0.22071175 0.21529464 - 0.21472976 + 0.21472977 0.20628634 0.2043854 - 0.20099047 + 0.20099048 0.18452412 0.17380331 0.15045828 - 0.13115592 + 0.13115594 0.1251017 0.10887553 0.079097189 0.061385516 0.043582208 - 0.026112026 + 0.026112027 0.01028904 -0.011082372 -0.031485982 @@ -216,43 +216,43 @@ -0.11735611 -0.12887773 -0.14123638 - -0.1465044 + -0.14650442 -0.16701357 -0.16918446 -0.17741001 - -0.18355082 + -0.18355083 -0.17439707 -0.17841232 -0.16846335 -0.16889346 - -0.17227842 + -0.17227843 -0.17070685 -0.15092036 -0.15334848 -0.12954405 -0.11295433 - -0.09564729 + -0.095647298 -0.080270432 -0.06479656 - -0.055495787 + -0.055495791 -0.03828676 -0.026680976 -0.0092150476 -0.0018779703 - 0.014258126 + 0.014258127 0.023644734 - 0.040185757 + 0.040185761 0.050761189 0.060287103 0.065608293 - 0.080362737 + 0.080362745 0.083356142 0.10414003 0.11332227 0.12000805 0.12286837 - 0.13463314 - 0.14336015 + 0.13463315 + 0.14336017 0.15017572 1.6947948e-05 1.720169e-05 @@ -505,5 +505,5 @@ -0.00016269117 -0.00016905692 - 0.049695708357404555 + 0.049695614293824519 diff --git a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP0.xml b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP0.xml index d318d513c0..6a528649ed 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP0.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP0.xml @@ -13,7 +13,7 @@ 0.65127754 0.58396322 0.50845456 - 0.43051061 + 0.43051064 0.35105249 0.26587203 0.18143819 @@ -43,7 +43,7 @@ -0.56959647 -0.53186774 -0.47233477 - -0.42289308 + -0.42289311 -0.36013916 -0.30238953 -0.23851043 @@ -116,7 +116,7 @@ 0.033095215 0.064681761 0.10605592 - 0.13676096 + 0.13676098 0.17619975 0.20634755 0.23424256 @@ -124,7 +124,7 @@ 0.28315401 0.31159228 0.32296148 - 0.3337082 + 0.33370823 0.3472271 0.35520938 0.35884556 @@ -167,7 +167,7 @@ -0.24391536 -0.23143035 -0.21355876 - -0.19516303 + -0.19516304 -0.17305651 -0.14851433 -0.13391782 @@ -178,7 +178,7 @@ -0.0082916291 0.028239559 0.044842608 - 0.066723041 + 0.066723049 0.086818486 0.11206019 0.13034853 @@ -201,7 +201,7 @@ 0.16507152 0.14102207 0.13159586 - 0.11218439 + 0.1121844 0.090035476 0.069503993 0.054288179 @@ -210,7 +210,7 @@ -0.0047988398 -0.024740862 -0.042112309 - -0.057106286 + -0.05710629 -0.079133794 -0.08531253 -0.10183473 @@ -505,5 +505,5 @@ -8.4526189e-05 -0.00019142308 - 0.54968381075777017 + 0.54968385173596346 diff --git a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP1.xml b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP1.xml index f02426b8fc..3784c42df8 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP1.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP1.xml @@ -38,7 +38,7 @@ -0.47069913 -0.44109708 -0.44801748 - -0.44286934 + -0.44286937 -0.40633747 -0.38159758 -0.38471025 @@ -46,7 +46,7 @@ -0.2954666 -0.22040513 -0.20538236 - -0.17562716 + -0.17562717 -0.12331922 -0.10807745 -0.050659552 @@ -125,7 +125,7 @@ 0.33220929 0.33306259 0.33927098 - 0.37198529 + 0.37198532 0.417128 0.4154191 0.39847392 @@ -238,7 +238,7 @@ -0.13059393 -0.10386769 -0.075353824 - -0.065707266 + -0.065707274 -0.018131474 -0.014584285 0.0061852229 @@ -505,5 +505,5 @@ -0.00067634159 -0.0014037765 - -0.13156737015469844 + -0.13156739250644023 diff --git a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP2.xml b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP2.xml index e61b5e32d2..7a673d6cb7 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP2.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP2.xml @@ -6,187 +6,187 @@ 1 1.0026263 1.002813 - 1.0028353 - 1.00298 + 1.0028352 + 1.0029799 1.0031264 - 1.0033009 + 1.0033008 1.0034248 1.0034949 - 1.0036868 + 1.0036867 1.0038811 - 1.0039461 + 1.0039459 1.0039393 - 1.0039014 - 1.0041367 - 1.0040401 + 1.0039012 + 1.0041366 + 1.00404 1.0039788 1.0041522 1.004038 - 1.00398 - 1.0038664 - 1.0038596 - 1.00374 - 1.0035481 - 1.0035183 + 1.0039799 + 1.0038663 + 1.0038595 + 1.0037398 + 1.003548 + 1.0035182 1.003296 1.0031067 - 1.0031048 - 1.0029247 + 1.0031047 + 1.0029246 1.0029281 - 1.0027844 + 1.0027843 1.0024828 1.0026224 - 1.0028158 + 1.0028157 1.0028605 1.0029123 - 1.0031923 + 1.0031922 1.0032749 - 1.0033575 - 1.0034987 - 1.0035222 - 1.0038122 + 1.0033574 + 1.0034986 + 1.003522 + 1.0038121 1.0038934 1.0038835 - 1.003985 + 1.0039849 1.0040591 1.0040226 1.004043 - 1.0039675 - 1.0040811 + 1.0039674 + 1.004081 1.0039179 - 1.003891 + 1.0038909 1.0038245 1.0036744 1.0038064 1.0035872 1.003276 - 1.0033685 + 1.0033684 1.0032185 1.0030403 1.0028248 1.0028371 - 1.0028143 + 1.0028142 1.00263 - 1.002896 + 1.0028958 1.0029092 - 1.0029544 + 1.0029542 1.0030841 1.0033208 - 1.0033398 - 1.0034574 - 1.003602 - 1.0037988 + 1.0033396 + 1.0034573 + 1.0036019 + 1.0037987 1.0038297 1.0040625 1.0039389 1.0040179 1.0039594 1.0039961 - 1.0039948 + 1.0039947 1.0040424 - 1.0039783 - 1.0039077 - 1.003881 - 1.003885 - 1.0038402 + 1.0039781 + 1.0039076 + 1.0038809 + 1.0038849 + 1.0038401 1.0036588 1.0035746 - 1.0034667 - 1.0033582 + 1.0034666 + 1.0033581 1.0032606 - 1.0031298 - 1.0029531 - 1.0029762 + 1.0031297 + 1.0029529 + 1.0029761 1.0028747 1.0028256 - 1.0029508 + 1.0029507 1.0031258 - 1.0032694 + 1.0032693 1.0033274 1.0034236 1.0035467 1.0037386 - 1.0037678 + 1.0037677 1.0038179 1.003939 1.003979 1.00396 1.0039968 - 1.0040765 + 1.0040764 1.0039867 1.0040717 1.0039583 1.0038927 - 1.0039638 + 1.0039637 1.0038794 1.003749 1.0036877 - 1.0035897 - 1.0034268 + 1.0035896 + 1.0034267 1.0034192 1.0032343 - 1.0033185 - 1.0029762 + 1.0033184 + 1.0029761 1.0031357 1.003207 - 1.0029137 + 1.0029136 1.0030683 - 1.0031815 + 1.0031813 1.0033453 - 1.0033225 + 1.0033224 1.0034435 - 1.0035906 - 1.0036092 + 1.0035905 + 1.0036091 1.0037911 - 1.0038091 - 1.0039287 + 1.003809 + 1.0039285 1.0040594 1.003963 1.0039438 - 1.004236 + 1.0042359 1.0041121 1.0039667 - 1.0039054 - 1.0039448 + 1.0039053 + 1.0039446 1.0038875 - 1.0039566 - 1.0037364 + 1.0039564 + 1.0037363 1.0038453 - 1.0036628 + 1.0036627 1.0036311 - 1.0035908 + 1.0035907 1.0035497 - 1.0035697 - 1.0033984 + 1.0035696 + 1.0033983 1.0034442 1.0033016 - 1.0032938 - 1.0032138 + 1.0032936 + 1.0032136 1.0033556 1.0033462 1.0034832 1.0035497 1.0036097 1.003722 - 1.003818 + 1.0038179 1.0037833 - 1.0038013 - 1.0039788 + 1.0038012 + 1.0039787 1.0039507 1.0040078 - 1.0040756 - 1.0040264 + 1.0040755 + 1.0040263 1.0038471 1.0039368 1.0040127 - 1.0038276 - 1.0038612 - 1.0038606 + 1.0038275 + 1.0038611 + 1.0038605 1.0038905 1.0039259 - 1.0036864 + 1.0036863 1.0036433 - 1.0036565 + 1.0036564 1.0036999 1.003649 1.0035149 @@ -195,21 +195,21 @@ 1.003474 1.0034895 1.003518 - 1.0036062 - 1.0036708 + 1.0036061 + 1.0036707 1.0038066 - 1.0037796 - 1.0038444 + 1.0037795 + 1.0038443 1.0038437 1.0038533 1.0040073 1.0039427 - 1.0040594 + 1.0040593 1.0040133 1.0039896 1.0041509 - 1.0040509 - 1.0040835 + 1.0040507 + 1.0040834 1.0039433 1.0039626 1.0039417 @@ -218,41 +218,41 @@ 1.0038795 1.0038284 1.0037864 - 1.0037988 - 1.0038092 - 1.0038675 - 1.0038058 + 1.0037987 + 1.0038091 + 1.0038674 + 1.0038056 1.0036768 - 1.0036308 + 1.0036306 1.0037299 - 1.003732 - 1.0038123 - 1.0038664 + 1.0037318 + 1.0038122 + 1.0038663 1.0037434 1.0037843 - 1.0037801 + 1.00378 1.0039681 - 1.003965 + 1.0039649 1.0039586 1.0040132 1.0040785 1.004027 1.0039933 - 1.0040841 + 1.004084 1.0040706 - 1.0040671 - 1.0040392 - 1.003906 - 1.0039387 - 1.0039802 + 1.0040669 + 1.004039 + 1.0039059 + 1.0039386 + 1.00398 1.0039933 1.0039399 1.0039351 - 1.0039452 + 1.0039451 1.0038861 - 1.0039293 + 1.0039291 1.0039737 - 1.0039887 + 1.0039886 1.0038258 -0.49893174 -0.49886957 @@ -505,5 +505,5 @@ -0.49923372 -0.49851903 - 127.17098122835159 + 127.17096823453903 diff --git a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP4.xml b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP4.xml index d318d513c0..6a528649ed 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP4.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacP4.xml @@ -13,7 +13,7 @@ 0.65127754 0.58396322 0.50845456 - 0.43051061 + 0.43051064 0.35105249 0.26587203 0.18143819 @@ -43,7 +43,7 @@ -0.56959647 -0.53186774 -0.47233477 - -0.42289308 + -0.42289311 -0.36013916 -0.30238953 -0.23851043 @@ -116,7 +116,7 @@ 0.033095215 0.064681761 0.10605592 - 0.13676096 + 0.13676098 0.17619975 0.20634755 0.23424256 @@ -124,7 +124,7 @@ 0.28315401 0.31159228 0.32296148 - 0.3337082 + 0.33370823 0.3472271 0.35520938 0.35884556 @@ -167,7 +167,7 @@ -0.24391536 -0.23143035 -0.21355876 - -0.19516303 + -0.19516304 -0.17305651 -0.14851433 -0.13391782 @@ -178,7 +178,7 @@ -0.0082916291 0.028239559 0.044842608 - 0.066723041 + 0.066723049 0.086818486 0.11206019 0.13034853 @@ -201,7 +201,7 @@ 0.16507152 0.14102207 0.13159586 - 0.11218439 + 0.1121844 0.090035476 0.069503993 0.054288179 @@ -210,7 +210,7 @@ -0.0047988398 -0.024740862 -0.042112309 - -0.057106286 + -0.05710629 -0.079133794 -0.08531253 -0.10183473 @@ -505,5 +505,5 @@ -8.4526189e-05 -0.00019142308 - 0.54968381075777017 + 0.54968385173596346 diff --git a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacVector.xml b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacVector.xml index d318d513c0..6a528649ed 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacVector.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/AutocorrTest_EacVector.xml @@ -13,7 +13,7 @@ 0.65127754 0.58396322 0.50845456 - 0.43051061 + 0.43051064 0.35105249 0.26587203 0.18143819 @@ -43,7 +43,7 @@ -0.56959647 -0.53186774 -0.47233477 - -0.42289308 + -0.42289311 -0.36013916 -0.30238953 -0.23851043 @@ -116,7 +116,7 @@ 0.033095215 0.064681761 0.10605592 - 0.13676096 + 0.13676098 0.17619975 0.20634755 0.23424256 @@ -124,7 +124,7 @@ 0.28315401 0.31159228 0.32296148 - 0.3337082 + 0.33370823 0.3472271 0.35520938 0.35884556 @@ -167,7 +167,7 @@ -0.24391536 -0.23143035 -0.21355876 - -0.19516303 + -0.19516304 -0.17305651 -0.14851433 -0.13391782 @@ -178,7 +178,7 @@ -0.0082916291 0.028239559 0.044842608 - 0.066723041 + 0.066723049 0.086818486 0.11206019 0.13034853 @@ -201,7 +201,7 @@ 0.16507152 0.14102207 0.13159586 - 0.11218439 + 0.1121844 0.090035476 0.069503993 0.054288179 @@ -210,7 +210,7 @@ -0.0047988398 -0.024740862 -0.042112309 - -0.057106286 + -0.05710629 -0.079133794 -0.08531253 -0.10183473 @@ -505,5 +505,5 @@ -8.4526189e-05 -0.00019142308 - 0.54968381075777017 + 0.54968385173596346 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnERF.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnERF.xml index f968c0c25e..b9f7053a5e 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnERF.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnERF.xml @@ -3,9 +3,9 @@ 4 - 7.8450705641969325 - 2.6181457557146546 - -124.26093946666471 - 10.991394787222672 + 93.478249768350366 + 106.54191922850268 + 180.89848007846916 + 4.0862253105513728 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnERREST.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnERREST.xml index 05f9939df8..5247edc079 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnERREST.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnERREST.xml @@ -3,8 +3,8 @@ 3 - 3.7998256591705597 - 0.82432886141340755 - 1.5224332184919882 + 0.50082180382778319 + 1 + 103.32999340041971 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP1.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP1.xml index b39c37344a..fc4b9df324 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP1.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP1.xml @@ -3,6 +3,6 @@ 1 - 28.790295709638542 + 28.633791718066412 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP2.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP2.xml index e30564ea47..23ffce35c3 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP2.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP2.xml @@ -3,7 +3,7 @@ 2 - 38.159752015476776 - 0.7956683286568571 + 36.978317574617037 + 0.80802279377014719 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP5.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP5.xml index 9dd9a8bc4a..e56125557b 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP5.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP5.xml @@ -3,10 +3,10 @@ 5 - 0.4874505859136089 - 5.593512789216291 - 0.5859548866294173 - 50.639858390563305 - -0.002429583757534859 + 0.59327596683805051 + 16.243870118801162 + 0.34166095324008267 + 96.46496015267067 + 0.014656658285601946 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP7.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP7.xml index 20c587183b..fdf0e31a98 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP7.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP7.xml @@ -3,12 +3,12 @@ 7 - 0.5512395081799536 - 6.5183422896301373 - -0.67296121974350642 - 31.927951035981206 - 1.186708933085455 - 42.008654192804514 - -0.0015455112108527936 + 0.22969109737961602 + 1.6938073214026785 + 0.58047615485173931 + 23.514547265682463 + 0.25468630875040837 + 124.69588927433887 + 0.0088943612374425204 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP9.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP9.xml index e9ecccc216..9998eae0b8 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP9.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP9.xml @@ -3,14 +3,14 @@ 9 - 15.957817686967116 - 221006.78653735726 - 114.88630569000553 - -33734.142061366547 - 7.6914607460936715 - 1876.2303193471271 - 0.75265207702909753 - 24.572552704802391 - -138.40801944191745 + 0.22969086672835151 + 1.6938036292509138 + 0.44194868847867191 + 23.513035628000434 + 0.13852713430709759 + 23.519216989391605 + 0.25468685075723652 + 124.69561162954194 + 0.008894417740064841 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP3.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXPEXP.xml similarity index 63% rename from src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP3.xml rename to src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXPEXP.xml index 74d4d3acb1..942d1ba996 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXP3.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnEXPEXP.xml @@ -3,8 +3,8 @@ 3 - 49.669368440495347 - 0.58492964998546348 - 6.4814246273214478 + 6.7235795110785457 + 0.42329531456542407 + 50.286846364014501 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnPRES.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnPRES.xml index fe3d48e213..64ba55b60c 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnPRES.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnPRES.xml @@ -3,11 +3,11 @@ 6 - 0.2214172263311579 - 10.066519430470473 - 3.7707829256907712 - 0.72109681793857272 - 0.49886880345093437 - 1.000154618776814 + 0.49146807677176813 + 10.122470863238238 + 7.3785890028734542 + 1.2022517448320804 + 0.41325325292628245 + 5.8277122335091622 diff --git a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnVAC.xml b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnVAC.xml index 35f0f01323..5fe01975e1 100644 --- a/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnVAC.xml +++ b/src/gromacs/correlationfunctions/tests/refdata/ExpfitTest_EffnVAC.xml @@ -3,7 +3,7 @@ 2 - 1.3759034165632673 - 0.20054113811854915 + 0.61188112533079864 + 0.085941133945086262 diff --git a/src/gromacs/correlationfunctions/tests/testERF.xvg b/src/gromacs/correlationfunctions/tests/testERF.xvg new file mode 100644 index 0000000000..79b426f436 --- /dev/null +++ b/src/gromacs/correlationfunctions/tests/testERF.xvg @@ -0,0 +1,401 @@ +0 93.486577 +1 92.437952 +2 92.676899 +3 93.287708 +4 92.524602 +5 93.33982 +6 94.322249 +7 94.152165 +8 93.268429 +9 92.985186 +10 92.447793 +11 93.529461 +12 93.811387 +13 92.128049 +14 93.058745 +15 91.834633 +16 94.592703 +17 94.835931 +18 92.57085 +19 93.350492 +20 93.793811 +21 93.064425 +22 92.358802 +23 93.138607 +24 91.978397 +25 92.883269 +26 94.889234 +27 94.240014 +28 94.141442 +29 94.588848 +30 93.267882 +31 92.642117 +32 94.11224 +33 93.724048 +34 93.050412 +35 92.166011 +36 94.357401 +37 92.104766 +38 93.609402 +39 94.193545 +40 92.616803 +41 93.190927 +42 94.279049 +43 94.066138 +44 92.766857 +45 93.945376 +46 93.410486 +47 92.710437 +48 94.479958 +49 93.37127 +50 94.423499 +51 94.57174 +52 93.645116 +53 93.276059 +54 94.777955 +55 92.01432 +56 93.443601 +57 94.161402 +58 92.767022 +59 94.703138 +60 92.866104 +61 92.709905 +62 93.443922 +63 92.749122 +64 92.996164 +65 93.226291 +66 94.187659 +67 92.239095 +68 93.12234 +69 93.177441 +70 92.728596 +71 92.900251 +72 94.832081 +73 92.066326 +74 93.740011 +75 93.962937 +76 93.559431 +77 92.615057 +78 94.143638 +79 94.78544 +80 93.031341 +81 94.033655 +82 92.418896 +83 92.221132 +84 92.18168 +85 93.919633 +86 94.881543 +87 93.927047 +88 94.674836 +89 93.190343 +90 92.947857 +91 93.298692 +92 94.621239 +93 94.459106 +94 92.906101 +95 94.237764 +96 93.645132 +97 93.660485 +98 92.915948 +99 93.317682 +100 92.120727 +101 93.230415 +102 92.674665 +103 92.159524 +104 93.868182 +105 92.996882 +106 94.363097 +107 94.502451 +108 92.358538 +109 93.717387 +110 92.786443 +111 92.186164 +112 92.229524 +113 92.36136 +114 94.888301 +115 94.425912 +116 92.947455 +117 92.82134 +118 94.760495 +119 94.977061 +120 93.748171 +121 93.605504 +122 94.126006 +123 94.857354 +124 94.460468 +125 92.21732 +126 93.804279 +127 93.52855 +128 94.313085 +129 94.035432 +130 92.303449 +131 94.071923 +132 92.482691 +133 95.028053 +134 92.937195 +135 92.615591 +136 93.580808 +137 94.951146 +138 92.656969 +139 93.256817 +140 92.455751 +141 95.206838 +142 92.984802 +143 92.945197 +144 94.699222 +145 94.56755 +146 95.347739 +147 92.459357 +148 93.171382 +149 95.180269 +150 94.758702 +151 93.379194 +152 93.86194 +153 92.852474 +154 93.279239 +155 93.569225 +156 94.28658 +157 94.625002 +158 93.888688 +159 94.990531 +160 94.31779 +161 93.264266 +162 95.189436 +163 95.923123 +164 93.687356 +165 94.856297 +166 94.09723 +167 96.538322 +168 95.588745 +169 94.174152 +170 94.815074 +171 96.021967 +172 95.894722 +173 97.48244 +174 95.421422 +175 95.738872 +176 96.590132 +177 98.483242 +178 98.857992 +179 97.832445 +180 99.176752 +181 102.0786 +182 101.71236 +183 100.84099 +184 102.06052 +185 103.04261 +186 104.13106 +187 102.48869 +188 103.97843 +189 103.45589 +190 103.46268 +191 103.76516 +192 104.77479 +193 104.83857 +194 103.92016 +195 105.16809 +196 104.87 +197 104.73319 +198 104.75198 +199 106.9468 +200 104.60863 +201 104.743 +202 106.03671 +203 104.92755 +204 104.56418 +205 105.0332 +206 106.65869 +207 104.39706 +208 106.46331 +209 104.91818 +210 106.79197 +211 105.90849 +212 105.78589 +213 106.76522 +214 106.13517 +215 106.84385 +216 107.33566 +217 105.98094 +218 105.46912 +219 104.71926 +220 106.89562 +221 106.95923 +222 106.8085 +223 104.86654 +224 105.4405 +225 106.28891 +226 107.53511 +227 106.43512 +228 105.04485 +229 107.42051 +230 104.91529 +231 105.72297 +232 106.4432 +233 104.84922 +234 107.01712 +235 107.34746 +236 107.16143 +237 106.14437 +238 105.83518 +239 105.85725 +240 106.36342 +241 107.85743 +242 107.0486 +243 105.25835 +244 106.05892 +245 105.21233 +246 105.15961 +247 104.99114 +248 104.98744 +249 106.47118 +250 107.77846 +251 105.95349 +252 106.0645 +253 105.00274 +254 106.7123 +255 106.78758 +256 107.71422 +257 106.74081 +258 107.91838 +259 107.14719 +260 107.57178 +261 107.71677 +262 106.37161 +263 107.76647 +264 105.35725 +265 107.10101 +266 106.95543 +267 106.38286 +268 107.87643 +269 106.18711 +270 105.68704 +271 107.61337 +272 107.30174 +273 107.55044 +274 107.66251 +275 107.97707 +276 108.02862 +277 105.89034 +278 105.8245 +279 106.14353 +280 106.56664 +281 105.80571 +282 107.29199 +283 105.60915 +284 106.63587 +285 106.82268 +286 107.21355 +287 106.76602 +288 106.24447 +289 106.10919 +290 105.73869 +291 105.88343 +292 107.37134 +293 105.92277 +294 107.80709 +295 106.73378 +296 105.28703 +297 107.99916 +298 107.20555 +299 106.72416 +300 106.3867 +301 106.73274 +302 106.99679 +303 107.56717 +304 106.47182 +305 107.21287 +306 106.39991 +307 107.90715 +308 105.09339 +309 107.98186 +310 107.65834 +311 107.39825 +312 107.62341 +313 106.49094 +314 105.92954 +315 107.08428 +316 105.46424 +317 106.1531 +318 105.92754 +319 107.00333 +320 107.87421 +321 105.98247 +322 105.65578 +323 106.11922 +324 105.21628 +325 105.15511 +326 107.23247 +327 107.6705 +328 107.93251 +329 105.8164 +330 105.19615 +331 107.22422 +332 106.46666 +333 105.1668 +334 107.92811 +335 106.38702 +336 106.0796 +337 106.50508 +338 106.15072 +339 105.26819 +340 105.40705 +341 107.25244 +342 105.59738 +343 106.33014 +344 106.44472 +345 106.88026 +346 107.05543 +347 105.61685 +348 106.21632 +349 106.35907 +350 107.49758 +351 106.7012 +352 106.10456 +353 106.93137 +354 107.94109 +355 106.36516 +356 105.22822 +357 107.60685 +358 107.94253 +359 107.26072 +360 105.4247 +361 106.41867 +362 107.69953 +363 105.63072 +364 105.9897 +365 106.46632 +366 106.01961 +367 106.09736 +368 107.5018 +369 106.07223 +370 107.5276 +371 107.19568 +372 107.08564 +373 106.88366 +374 105.74583 +375 107.09943 +376 107.08013 +377 107.45954 +378 107.53612 +379 108.112 +380 107.86845 +381 107.9672 +382 105.9697 +383 107.29257 +384 107.10105 +385 105.21125 +386 107.27416 +387 107.42193 +388 106.61502 +389 106.68435 +390 106.92447 +391 106.36125 +392 106.20247 +393 106.93532 +394 106.22898 +395 106.99196 +396 106.53996 +397 107.53919 +398 105.43453 +399 107.28519 +400 105.87991 diff --git a/src/gromacs/correlationfunctions/tests/testERREST.xvg b/src/gromacs/correlationfunctions/tests/testERREST.xvg new file mode 100644 index 0000000000..4f28dfc34b --- /dev/null +++ b/src/gromacs/correlationfunctions/tests/testERREST.xvg @@ -0,0 +1,501 @@ +1 0.3517423 +17 1.0597503 +33 0.9410293 +49 0.8716023 +65 0.8238423 +81 1.0098463 +97 0.5520603 +113 0.4278303 +129 0.9342923 +145 1.0342603 +161 0.8231863 +177 0.7635853 +193 0.9390683 +209 0.8922163 +225 1.1775333 +241 1.2116213 +257 0.6996663 +273 0.2287793 +289 1.3589063 +305 1.6198323 +321 0.8974273 +337 0.9043013 +353 0.4556703 +369 0.6929363 +385 1.2148863 +401 1.4234283 +417 1.0278523 +433 0.6124313 +449 0.7412923 +465 0.3195843 +481 0.8295103 +497 0.7101183 +513 1.0936703 +529 0.4371843 +545 0.5090993 +561 0.9473693 +577 0.9137463 +593 0.6570103 +609 1.0386173 +625 1.3944973 +641 0.4820373 +657 0.4620413 +673 1.3725933 +689 0.8170593 +705 0.8616223 +721 0.4220783 +737 1.1177493 +753 0.7920653 +769 0.6135533 +785 0.6946533 +801 1.4141893 +817 0.9795273 +833 0.5320023 +849 1.1125843 +865 1.2062423 +881 0.9650543 +897 1.6114013 +913 0.5465363 +929 1.5686003 +945 1.3882573 +961 1.0649153 +977 1.3672613 +993 1.1256913 +1009 0.8110243 +1025 1.0385253 +1041 1.4888273 +1057 1.4508483 +1073 0.9516873 +1089 0.9449663 +1105 0.9077883 +1121 0.5280203 +1137 0.8003353 +1153 0.3715183 +1169 1.1737183 +1185 0.3641943 +1201 0.9924893 +1217 1.4245123 +1233 0.7253083 +1249 1.0137533 +1265 0.8326003 +1281 0.8174633 +1297 0.6439253 +1313 0.9496813 +1329 0.4628573 +1345 1.3082933 +1361 0.3596463 +1377 0.8700153 +1393 1.0316053 +1409 1.2754943 +1425 0.8401303 +1441 1.3051733 +1457 1.7206623 +1473 0.9644593 +1489 0.7674533 +1505 0.9751783 +1521 0.5897953 +1537 0.9886143 +1553 1.5970503 +1569 1.3752943 +1585 1.1988803 +1601 1.0081603 +1617 1.5943493 +1633 0.7405973 +1649 0.6524783 +1665 0.7423293 +1681 1.1931963 +1697 1.2285283 +1713 0.7597933 +1729 0.9291963 +1745 1.5914353 +1761 0.4776883 +1777 1.3845643 +1793 1.4108173 +1809 0.8124203 +1825 1.4891713 +1841 0.8480123 +1857 0.7482803 +1873 1.5186353 +1889 0.8382693 +1905 0.3961233 +1921 0.8167013 +1937 1.3302663 +1953 0.5777863 +1969 1.4406173 +1985 0.4306763 +2001 0.6667753 +2017 -0.0174287 +2033 1.0792663 +2049 0.8624543 +2065 0.7621893 +2081 0.7339673 +2097 0.5926183 +2113 1.0041013 +2129 0.5759933 +2145 1.4805423 +2161 0.9656263 +2177 0.7245373 +2193 0.7458013 +2209 0.7898533 +2225 1.5943193 +2241 0.6717653 +2257 1.2236373 +2273 0.5686843 +2289 1.4059033 +2305 0.8901873 +2321 1.1125003 +2337 0.9947173 +2353 0.7398113 +2369 1.4866303 +2385 1.7170683 +2401 1.1582763 +2417 0.4476063 +2433 0.8380023 +2449 0.7195173 +2465 1.4064993 +2481 1.1095623 +2497 8.13e-05 +2513 0.8899893 +2529 0.9137923 +2545 1.0530973 +2561 1.2563603 +2577 0.3981443 +2593 0.8871053 +2609 0.4778103 +2625 1.1359983 +2641 0.5740863 +2657 1.1778153 +2673 1.3080033 +2689 1.0777173 +2705 1.7451983 +2721 1.9104733 +2737 1.4288603 +2753 2.2134143 +2769 0.9296233 +2785 0.9103213 +2801 0.8418393 +2817 1.2955373 +2833 1.2594723 +2849 0.8617293 +2865 0.2954833 +2881 0.9188513 +2897 1.8631033 +2913 1.9094973 +2929 0.7062733 +2945 1.8554883 +2961 0.9238253 +2977 0.3348663 +2993 0.9353073 +3009 1.2028623 +3025 1.9851123 +3041 1.1466263 +3057 1.5713393 +3073 1.4108013 +3089 0.8622563 +3105 0.8434723 +3121 1.3925823 +3137 1.4874463 +3153 1.2751583 +3169 1.1509373 +3185 1.5736133 +3201 1.5143093 +3217 0.7679413 +3233 1.0632833 +3249 0.9187743 +3265 1.5776033 +3281 1.1580243 +3297 0.8602263 +3313 0.8002063 +3329 0.9148603 +3345 1.3152133 +3361 0.4933133 +3377 0.8761343 +3393 0.6942183 +3409 1.2085083 +3425 1.5972033 +3441 1.2005583 +3457 1.2183353 +3473 0.7577103 +3489 0.6112263 +3505 1.2318773 +3521 1.1326723 +3537 1.1506243 +3553 0.5817003 +3569 0.9014713 +3585 0.8777513 +3601 1.4347963 +3617 0.9277923 +3633 1.1432003 +3649 0.7949343 +3665 1.2829793 +3681 0.2914853 +3697 0.9427083 +3713 1.3510713 +3729 0.8358273 +3745 0.8312423 +3761 1.0880323 +3777 1.7421923 +3793 1.1334733 +3809 0.2839173 +3825 0.8868153 +3841 0.5585983 +3857 0.5283023 +3873 0.2076003 +3889 0.8009313 +3905 0.9836093 +3921 0.8874943 +3937 0.5242053 +3953 0.6777083 +3969 0.2986573 +3985 0.9648863 +4001 0.3936813 +4017 1.5192383 +4033 1.1283383 +4049 0.5402493 +4065 0.8644763 +4081 0.9670683 +4097 1.2198303 +4113 0.8354763 +4129 0.8359803 +4145 1.0664183 +4161 0.4950833 +4177 0.6145223 +4193 0.8997543 +4209 0.2690933 +4225 0.9369703 +4241 1.3704423 +4257 0.5915883 +4273 0.7297033 +4289 0.2907223 +4305 1.0001723 +4321 1.1066633 +4337 1.5133253 +4353 1.3272213 +4369 1.2639593 +4385 0.6724293 +4401 0.6174593 +4417 1.4950913 +4433 0.9158143 +4449 1.9446683 +4465 1.1331683 +4481 0.5225033 +4497 0.8273663 +4513 0.3377963 +4529 0.6311003 +4545 1.4785433 +4561 0.8769733 +4577 0.7451523 +4593 1.0181933 +4609 1.1503113 +4625 1.5157673 +4641 1.2449163 +4657 1.5840353 +4673 1.1181763 +4689 0.8027923 +4705 0.2972533 +4721 0.6047333 +4737 1.3964663 +4753 0.8300983 +4769 1.5832113 +4785 1.7844893 +4801 1.4584933 +4817 1.2985273 +4833 1.1074343 +4849 1.6353573 +4865 1.5705303 +4881 1.7136353 +4897 1.3027233 +4913 1.0545013 +4929 1.0732243 +4945 1.9772463 +4961 1.4843033 +4977 1.4208723 +4993 1.8630953 +5009 0.8611723 +5025 0.7367903 +5041 0.5593233 +5057 0.9989823 +5073 0.5387463 +5089 0.6843843 +5105 0.7769823 +5121 1.1148653 +5137 1.7723363 +5153 0.6389133 +5169 0.6624193 +5185 0.5881323 +5201 0.5588963 +5217 -0.1075927 +5233 0.5866443 +5249 0.4331333 +5265 0.7912493 +5281 0.1760223 +5297 0.8659643 +5313 1.1929063 +5329 0.4952133 +5345 0.6245163 +5361 0.6433383 +5377 0.6050383 +5393 1.1802343 +5409 0.9662523 +5425 1.1449483 +5441 0.4728213 +5457 0.7832533 +5473 1.0352293 +5489 0.8087053 +5505 0.6213123 +5521 0.8185243 +5537 0.2622273 +5553 0.8098043 +5569 0.3489813 +5585 0.5309423 +5601 0.5585833 +5617 -0.0650127 +5633 0.8007403 +5649 1.1783113 +5665 0.8767213 +5681 0.7084473 +5697 0.7087373 +5713 1.2293673 +5729 1.9700823 +5745 1.1370283 +5761 0.8622563 +5777 1.5045673 +5793 1.5535933 +5809 1.6912733 +5825 1.4561283 +5841 1.1133853 +5857 1.0418363 +5873 1.4926883 +5889 0.7510653 +5905 1.4871873 +5921 0.6951643 +5937 0.8722043 +5953 1.0236633 +5969 0.4982503 +5985 1.0484053 +6001 1.4949693 +6017 0.7577103 +6033 1.4756743 +6049 0.9504673 +6065 0.6371503 +6081 0.9289373 +6097 1.6262783 +6113 1.5682493 +6129 1.7294363 +6145 0.8229643 +6161 1.4206823 +6177 1.2592053 +6193 0.8906453 +6209 0.9228333 +6225 0.7128493 +6241 1.0956393 +6257 1.3999223 +6273 1.2278183 +6289 1.3250013 +6305 0.7005813 +6321 0.8609823 +6337 0.9442873 +6353 0.7324343 +6369 1.2111333 +6385 0.9367413 +6401 0.5429203 +6417 0.7491193 +6433 0.6651583 +6449 1.0823103 +6465 1.3247193 +6481 0.8797353 +6497 1.0703243 +6513 0.9107563 +6529 1.8092243 +6545 0.7302523 +6561 1.6127293 +6577 1.2735493 +6593 1.0610703 +6609 0.7459463 +6625 1.0303013 +6641 1.5765423 +6657 1.0391893 +6673 1.1801883 +6689 0.3454793 +6705 0.5029193 +6721 0.7490663 +6737 1.1911363 +6753 1.1398823 +6769 0.8100933 +6785 1.4742703 +6801 1.1287503 +6817 1.4290363 +6833 1.2512253 +6849 1.1904803 +6865 1.1241423 +6881 0.9222003 +6897 0.5952273 +6913 0.8656583 +6929 1.1323363 +6945 1.1635483 +6961 1.3734173 +6977 1.2152753 +6993 0.7563443 +7009 0.6301693 +7025 0.3757603 +7041 0.9500623 +7057 0.9807333 +7073 0.6007513 +7089 1.1412473 +7105 0.8297853 +7121 0.8890813 +7137 0.8713803 +7153 1.0812503 +7169 1.3144803 +7185 1.3724333 +7201 1.7222793 +7217 0.8696183 +7233 1.2517673 +7249 1.0543183 +7265 1.0174603 +7281 1.1306963 +7297 0.4270143 +7313 0.9045993 +7329 1.5317583 +7345 1.1389813 +7361 0.6676073 +7377 0.9130753 +7393 0.5875443 +7409 0.9053693 +7425 0.5835163 +7441 1.1831793 +7457 1.4384203 +7473 0.4676183 +7489 0.7656293 +7505 0.4548773 +7521 1.0648013 +7537 1.0835923 +7553 1.3808793 +7569 1.2125673 +7585 1.2858243 +7601 0.7112093 +7617 0.6617253 +7633 0.8534213 +7649 1.2384923 +7665 1.1525693 +7681 1.8521543 +7697 1.2422913 +7713 1.4445463 +7729 1.0739483 +7745 0.8565343 +7761 0.6647313 +7777 1.2838183 +7793 1.6710863 +7809 1.7629593 +7825 1.7022833 +7841 1.7739233 +7857 1.3768663 +7873 1.7267503 +7889 1.3419463 +7905 1.4587673 +7921 0.8886153 +7937 0.6786163 +7953 0.9597973 +7969 1.7429853 +7985 1.1981473 +8001 1.7303743 diff --git a/src/gromacs/correlationfunctions/tests/testEXP.xvg b/src/gromacs/correlationfunctions/tests/testEXP.xvg deleted file mode 100644 index a192c51c80..0000000000 --- a/src/gromacs/correlationfunctions/tests/testEXP.xvg +++ /dev/null @@ -1,501 +0,0 @@ -0 0 -1 0.933793 -2 1.64413 -3 2.26881 -4 2.69732 -5 3.18609 -6 3.5633 -7 3.77591 -8 4.083 -9 4.27165 -10 4.50143 -11 4.68664 -12 4.84612 -13 4.93256 -14 5.0487 -15 5.21257 -16 5.29916 -17 5.32257 -18 5.43037 -19 5.50601 -20 5.60239 -21 5.62398 -22 5.72622 -23 5.77308 -24 5.7562 -25 5.84888 -26 5.81699 -27 5.86178 -28 5.94235 -29 5.95981 -30 5.98369 -31 5.95674 -32 6.06792 -33 6.07991 -34 6.11024 -35 6.12377 -36 6.06567 -37 6.16848 -38 6.15328 -39 6.16295 -40 6.16469 -41 6.22697 -42 6.1776 -43 6.27148 -44 6.23865 -45 6.2744 -46 6.27211 -47 6.29498 -48 6.31175 -49 6.28439 -50 6.34758 -51 6.26659 -52 6.29551 -53 6.31016 -54 6.31924 -55 6.32594 -56 6.39353 -57 6.37555 -58 6.4125 -59 6.34882 -60 6.34253 -61 6.37401 -62 6.42251 -63 6.44891 -64 6.41029 -65 6.36925 -66 6.45488 -67 6.47461 -68 6.39323 -69 6.41703 -70 6.47664 -71 6.47864 -72 6.40653 -73 6.4444 -74 6.4901 -75 6.49016 -76 6.43924 -77 6.50607 -78 6.48499 -79 6.4465 -80 6.51673 -81 6.44586 -82 6.53225 -83 6.51432 -84 6.52803 -85 6.49069 -86 6.51123 -87 6.54134 -88 6.50065 -89 6.49403 -90 6.54445 -91 6.56578 -92 6.54692 -93 6.58151 -94 6.5848 -95 6.5796 -96 6.53327 -97 6.57402 -98 6.51995 -99 6.59609 -100 6.59545 -101 6.51956 -102 6.54083 -103 6.53558 -104 6.51182 -105 6.5347 -106 6.55378 -107 6.58284 -108 6.56471 -109 6.57831 -110 6.59092 -111 6.62654 -112 6.54307 -113 6.57412 -114 6.60729 -115 6.59664 -116 6.63216 -117 6.59422 -118 6.61469 -119 6.62636 -120 6.57808 -121 6.57356 -122 6.63797 -123 6.58006 -124 6.62347 -125 6.56521 -126 6.63422 -127 6.57666 -128 6.56846 -129 6.57003 -130 6.56552 -131 6.56682 -132 6.58905 -133 6.59197 -134 6.58329 -135 6.57486 -136 6.6287 -137 6.62418 -138 6.57194 -139 6.62585 -140 6.59216 -141 6.66186 -142 6.61182 -143 6.64493 -144 6.60235 -145 6.64946 -146 6.58789 -147 6.636 -148 6.60824 -149 6.59174 -150 6.67204 -151 6.63844 -152 6.64634 -153 6.61669 -154 6.60086 -155 6.68471 -156 6.6902 -157 6.65257 -158 6.61755 -159 6.67101 -160 6.6674 -161 6.63823 -162 6.63133 -163 6.6367 -164 6.66352 -165 6.696 -166 6.60881 -167 6.67253 -168 6.66735 -169 6.66583 -170 6.66259 -171 6.63691 -172 6.6766 -173 6.6913 -174 6.65969 -175 6.61821 -176 6.6796 -177 6.70715 -178 6.67694 -179 6.69653 -180 6.64676 -181 6.65935 -182 6.6781 -183 6.71181 -184 6.70878 -185 6.67739 -186 6.64651 -187 6.66021 -188 6.71047 -189 6.65987 -190 6.6511 -191 6.71258 -192 6.71689 -193 6.65441 -194 6.65054 -195 6.66806 -196 6.70738 -197 6.64787 -198 6.69967 -199 6.6331 -200 6.63513 -201 6.69851 -202 6.64237 -203 6.70702 -204 6.70229 -205 6.66568 -206 6.6818 -207 6.68629 -208 6.71863 -209 6.73129 -210 6.66952 -211 6.69274 -212 6.67642 -213 6.6782 -214 6.65372 -215 6.7269 -216 6.70297 -217 6.64193 -218 6.65786 -219 6.73524 -220 6.70002 -221 6.6906 -222 6.70867 -223 6.7153 -224 6.66604 -225 6.70376 -226 6.70701 -227 6.71955 -228 6.73206 -229 6.71019 -230 6.70016 -231 6.72956 -232 6.64624 -233 6.73931 -234 6.70839 -235 6.67605 -236 6.67944 -237 6.7271 -238 6.72234 -239 6.66538 -240 6.71936 -241 6.66755 -242 6.68937 -243 6.724 -244 6.68827 -245 6.69197 -246 6.66415 -247 6.69317 -248 6.69017 -249 6.70377 -250 6.70761 -251 6.74691 -252 6.7327 -253 6.66325 -254 6.67565 -255 6.68625 -256 6.71185 -257 6.65938 -258 6.66814 -259 6.69835 -260 6.71306 -261 6.68137 -262 6.68514 -263 6.69503 -264 6.70687 -265 6.70667 -266 6.71226 -267 6.74293 -268 6.6743 -269 6.70273 -270 6.68393 -271 6.71818 -272 6.72623 -273 6.73187 -274 6.73242 -275 6.70873 -276 6.70823 -277 6.69239 -278 6.75188 -279 6.71371 -280 6.7502 -281 6.69877 -282 6.6998 -283 6.66478 -284 6.70866 -285 6.73126 -286 6.73144 -287 6.70046 -288 6.70388 -289 6.6847 -290 6.75264 -291 6.73971 -292 6.70705 -293 6.7407 -294 6.67035 -295 6.72357 -296 6.73306 -297 6.73084 -298 6.68237 -299 6.76651 -300 6.72838 -301 6.76579 -302 6.744 -303 6.6963 -304 6.76216 -305 6.66966 -306 6.67184 -307 6.67226 -308 6.67163 -309 6.68679 -310 6.73726 -311 6.68492 -312 6.73491 -313 6.75372 -314 6.68562 -315 6.73127 -316 6.73742 -317 6.74772 -318 6.71254 -319 6.67357 -320 6.72034 -321 6.71205 -322 6.68542 -323 6.76568 -324 6.69664 -325 6.7009 -326 6.71605 -327 6.73142 -328 6.7524 -329 6.70307 -330 6.7459 -331 6.70831 -332 6.71814 -333 6.76659 -334 6.67721 -335 6.70581 -336 6.75936 -337 6.71626 -338 6.77537 -339 6.73916 -340 6.75086 -341 6.69825 -342 6.73784 -343 6.68206 -344 6.73098 -345 6.69878 -346 6.70428 -347 6.76935 -348 6.68665 -349 6.72539 -350 6.70694 -351 6.683 -352 6.76606 -353 6.74774 -354 6.69002 -355 6.77377 -356 6.71343 -357 6.71858 -358 6.77459 -359 6.72939 -360 6.74046 -361 6.71403 -362 6.70104 -363 6.7676 -364 6.70533 -365 6.68404 -366 6.73426 -367 6.74651 -368 6.73703 -369 6.7093 -370 6.68535 -371 6.78054 -372 6.68628 -373 6.74885 -374 6.75663 -375 6.78102 -376 6.78231 -377 6.69437 -378 6.69903 -379 6.70758 -380 6.74678 -381 6.7422 -382 6.72993 -383 6.68776 -384 6.74261 -385 6.75407 -386 6.76879 -387 6.70534 -388 6.6894 -389 6.71274 -390 6.73197 -391 6.75648 -392 6.74171 -393 6.73153 -394 6.68943 -395 6.75898 -396 6.72838 -397 6.74778 -398 6.72736 -399 6.77489 -400 6.77639 -401 6.70282 -402 6.71161 -403 6.76628 -404 6.77065 -405 6.77603 -406 6.68948 -407 6.76621 -408 6.72207 -409 6.72234 -410 6.75854 -411 6.77273 -412 6.76151 -413 6.74252 -414 6.73153 -415 6.69708 -416 6.74061 -417 6.75818 -418 6.72091 -419 6.74424 -420 6.78673 -421 6.74145 -422 6.69344 -423 6.72471 -424 6.72518 -425 6.76767 -426 6.71893 -427 6.78251 -428 6.78466 -429 6.78094 -430 6.77212 -431 6.74085 -432 6.74283 -433 6.73333 -434 6.78714 -435 6.77238 -436 6.73692 -437 6.73205 -438 6.70148 -439 6.77486 -440 6.76451 -441 6.79155 -442 6.73363 -443 6.69988 -444 6.77442 -445 6.78397 -446 6.7856 -447 6.73778 -448 6.70429 -449 6.7378 -450 6.74829 -451 6.69527 -452 6.75622 -453 6.7867 -454 6.76332 -455 6.72525 -456 6.72494 -457 6.78544 -458 6.7868 -459 6.73629 -460 6.75375 -461 6.71361 -462 6.73621 -463 6.74825 -464 6.70532 -465 6.78782 -466 6.73666 -467 6.79325 -468 6.69962 -469 6.71673 -470 6.72258 -471 6.72758 -472 6.70651 -473 6.70794 -474 6.77938 -475 6.71007 -476 6.7414 -477 6.77997 -478 6.71039 -479 6.77047 -480 6.71826 -481 6.73399 -482 6.69929 -483 6.7929 -484 6.77405 -485 6.79742 -486 6.76882 -487 6.74596 -488 6.71524 -489 6.7947 -490 6.75326 -491 6.72172 -492 6.71176 -493 6.71697 -494 6.74455 -495 6.78932 -496 6.70312 -497 6.72368 -498 6.70806 -499 6.70917 -500 6.79967 diff --git a/src/gromacs/correlationfunctions/tests/testINVEXP.xvg b/src/gromacs/correlationfunctions/tests/testINVEXP.xvg old mode 100755 new mode 100644 diff --git a/src/gromacs/correlationfunctions/tests/testINVEXP79.xvg b/src/gromacs/correlationfunctions/tests/testINVEXP79.xvg new file mode 100644 index 0000000000..f96379f7d5 --- /dev/null +++ b/src/gromacs/correlationfunctions/tests/testINVEXP79.xvg @@ -0,0 +1,501 @@ +0 1.09007 +1 0.877207 +2 0.932819 +3 0.810603 +4 0.791896 +5 0.733364 +6 0.616033 +7 0.604703 +8 0.755853 +9 0.636944 +10 0.715246 +11 0.611031 +12 0.595407 +13 0.63075 +14 0.581412 +15 0.480107 +16 0.51343 +17 0.456788 +18 0.44948 +19 0.435668 +20 0.403863 +21 0.489832 +22 0.495195 +23 0.422594 +24 0.480228 +25 0.482537 +26 0.444401 +27 0.384561 +28 0.34734 +29 0.383506 +30 0.301289 +31 0.382879 +32 0.292707 +33 0.349764 +34 0.302061 +35 0.407986 +36 0.434076 +37 0.335145 +38 0.260494 +39 0.355107 +40 0.25965 +41 0.296108 +42 0.251938 +43 0.255597 +44 0.277503 +45 0.396491 +46 0.28228 +47 0.27551 +48 0.270726 +49 0.292784 +50 0.128415 +51 0.132945 +52 0.23486 +53 0.241426 +54 0.156775 +55 0.2457 +56 0.174162 +57 0.349717 +58 0.186471 +59 0.300814 +60 0.15326 +61 0.231873 +62 0.163605 +63 0.197539 +64 0.260759 +65 0.231805 +66 0.233305 +67 0.135304 +68 0.281354 +69 0.234155 +70 0.176055 +71 0.198226 +72 0.213007 +73 0.225185 +74 0.274865 +75 0.247751 +76 0.105759 +77 0.185042 +78 0.0998593 +79 0.19495 +80 0.0787285 +81 0.31299 +82 0.116831 +83 0.0640667 +84 0.111865 +85 0.127155 +86 0.0863288 +87 0.0788565 +88 0.229011 +89 0.0803561 +90 0.0798449 +91 0.103291 +92 0.148852 +93 0.0735688 +94 0.159655 +95 0.108183 +96 0.0602938 +97 0.148602 +98 0.148881 +99 0.115299 +100 0.0788461 +101 0.181537 +102 0.178337 +103 0.130516 +104 0.0405722 +105 0.178499 +106 0.156721 +107 0.117595 +108 0.190161 +109 0.0353169 +110 0.0107063 +111 0.044126 +112 0.179688 +113 0.0826383 +114 0.161516 +115 0.211201 +116 -0.0253689 +117 0.198065 +118 0.0845166 +119 0.0983947 +120 0.0836944 +121 0.00686212 +122 0.138417 +123 -0.0243976 +124 0.0290124 +125 0.192813 +126 0.190617 +127 0.153298 +128 0.153291 +129 0.160775 +130 0.0652794 +131 0.122779 +132 0.240448 +133 0.072505 +134 0.140685 +135 0.0979626 +136 -0.033676 +137 0.0573126 +138 0.162429 +139 0.0502877 +140 0.176007 +141 0.059104 +142 0.119028 +143 0.0682911 +144 0.0443568 +145 0.0813819 +146 0.046495 +147 0.213421 +148 0.116277 +149 0.211761 +150 0.008448 +151 0.143644 +152 0.203298 +153 0.0675614 +154 0.0783106 +155 0.0764201 +156 0.133036 +157 0.044048 +158 0.0901729 +159 0.00259126 +160 0.088098 +161 0.0468678 +162 0.0648707 +163 0.242762 +164 0.0509743 +165 0.0502586 +166 -0.00215295 +167 0.0977259 +168 0.211509 +169 0.122585 +170 0.115359 +171 0.124501 +172 0.0934846 +173 0.0607657 +174 0.0455892 +175 0.0295404 +176 -0.0230954 +177 0.115156 +178 0.171233 +179 0.0706772 +180 0.0125189 +181 0.0264475 +182 -0.0298884 +183 0.119367 +184 -0.0348317 +185 0.0161417 +186 0.210381 +187 0.15479 +188 0.0625029 +189 0.112174 +190 0.0869022 +191 0.101378 +192 0.0456388 +193 0.0512337 +194 0.0453277 +195 0.141329 +196 0.0583165 +197 -0.100644 +198 -0.0195298 +199 0.0536575 +200 0.0711223 +201 0.0323407 +202 0.0755175 +203 0.151374 +204 0.117336 +205 -0.0224392 +206 0.128654 +207 -0.0316783 +208 -0.00303418 +209 0.15873 +210 -0.0463493 +211 0.027255 +212 0.125538 +213 0.165183 +214 -0.0190088 +215 -0.0792838 +216 0.0241941 +217 0.027236 +218 0.0460868 +219 0.116271 +220 0.00766963 +221 0.0391876 +222 0.011154 +223 0.117667 +224 -0.0332919 +225 0.0515272 +226 -0.0726914 +227 0.0734784 +228 0.00147696 +229 0.133906 +230 0.0234433 +231 0.0647977 +232 0.162547 +233 0.0904404 +234 0.14845 +235 0.00353783 +236 0.175999 +237 0.0489034 +238 0.0733482 +239 0.052908 +240 0.00921323 +241 0.165374 +242 -0.0401146 +243 -0.0183761 +244 0.0339211 +245 -0.0140982 +246 0.162167 +247 0.083416 +248 0.0430574 +249 -0.0532527 +250 0.147883 +251 0.0333986 +252 -0.00200228 +253 0.11477 +254 0.105691 +255 -0.0071386 +256 0.04273 +257 0.0495165 +258 0.0653199 +259 0.0540464 +260 -0.0184714 +261 -0.0314251 +262 -0.0630476 +263 -0.0235271 +264 0.0366917 +265 0.0977608 +266 0.0564583 +267 0.106372 +268 0.0429041 +269 -0.0523068 +270 0.00302942 +271 0.0707971 +272 0.134789 +273 -0.00167497 +274 0.0757873 +275 -0.0220373 +276 0.00421086 +277 0.0529898 +278 0.0759217 +279 0.0753902 +280 0.103891 +281 0.023528 +282 0.0209249 +283 0.117407 +284 -0.00610032 +285 -0.0319453 +286 0.0774823 +287 -0.0464274 +288 0.0475878 +289 -0.0570683 +290 -0.0604167 +291 -0.0104142 +292 0.0735168 +293 0.0105054 +294 0.138414 +295 -0.0477548 +296 -0.0100453 +297 0.0441746 +298 -0.0494035 +299 0.0100155 +300 -0.0115031 +301 -0.0476504 +302 -0.0486697 +303 0.000769082 +304 0.0833754 +305 0.0828655 +306 0.0790107 +307 -0.0217383 +308 -0.0369806 +309 0.0294119 +310 -0.0621266 +311 0.0289071 +312 -0.0395707 +313 0.0415419 +314 0.085237 +315 0.050411 +316 -0.0210736 +317 0.0998225 +318 0.109902 +319 -0.0361716 +320 -0.0101586 +321 -0.0223105 +322 0.00826686 +323 -0.0293934 +324 0.0987037 +325 0.0315682 +326 0.104734 +327 0.0861752 +328 0.0764943 +329 0.000824613 +330 -0.0170506 +331 0.00164694 +332 0.0325245 +333 0.0791064 +334 0.11446 +335 -0.0148082 +336 0.0455031 +337 0.0586395 +338 -0.108668 +339 0.0346456 +340 0.0293549 +341 0.0890408 +342 -0.0106896 +343 0.112917 +344 -0.0296479 +345 0.0522096 +346 0.208064 +347 -0.00224504 +348 0.0121225 +349 -0.105139 +350 -0.00490689 +351 0.0974383 +352 0.045095 +353 0.0802632 +354 -0.000141307 +355 0.00395148 +356 0.156947 +357 0.0607425 +358 0.0246873 +359 0.0152081 +360 0.0806166 +361 0.0525148 +362 0.0390263 +363 -0.082116 +364 -0.0435646 +365 0.0144118 +366 -0.056396 +367 0.0224241 +368 0.0249577 +369 0.032059 +370 0.0710599 +371 0.0617842 +372 0.0870343 +373 -0.0309283 +374 -0.0503928 +375 0.0654705 +376 -0.132297 +377 -0.013171 +378 0.15629 +379 -0.0195602 +380 0.0804519 +381 0.0743904 +382 -0.013602 +383 -0.0495775 +384 -0.0715722 +385 0.150116 +386 0.00559995 +387 0.0252303 +388 0.0791371 +389 -0.00344736 +390 0.051361 +391 -0.0283868 +392 0.00921561 +393 0.0701514 +394 0.110203 +395 -0.00277509 +396 0.155823 +397 0.00737347 +398 -0.00262632 +399 -0.020737 +400 0.0512275 +401 -0.011247 +402 -0.0145163 +403 0.00934317 +404 0.157135 +405 0.102475 +406 0.0686777 +407 0.00617163 +408 -0.024658 +409 0.0290312 +410 0.0888993 +411 0.0475782 +412 -0.0887449 +413 0.068523 +414 -0.029187 +415 0.0134721 +416 0.0268544 +417 0.11216 +418 -0.0280536 +419 -0.0487176 +420 -0.000343852 +421 -0.0204678 +422 0.0346145 +423 -0.0845521 +424 -0.0557614 +425 0.123208 +426 -0.0138119 +427 -0.0674846 +428 -0.0271054 +429 0.0822262 +430 0.0342709 +431 -0.0765876 +432 0.0319721 +433 -0.0739633 +434 0.0233444 +435 -0.00554542 +436 0.0904659 +437 -0.00953448 +438 -0.0605533 +439 -0.0713493 +440 0.017423 +441 0.067253 +442 0.009445 +443 0.00886866 +444 0.0647234 +445 0.0687361 +446 0.00146581 +447 0.0044238 +448 0.0517999 +449 -0.110766 +450 0.0330447 +451 -0.0273471 +452 0.00885882 +453 -0.0352997 +454 0.133558 +455 0.107036 +456 0.0951639 +457 0.0922061 +458 0.107501 +459 0.00342741 +460 0.0689463 +461 -0.109189 +462 0.0251697 +463 -0.104037 +464 0.00550332 +465 0.0253927 +466 0.0473155 +467 -0.0858022 +468 0.0200478 +469 0.0549883 +470 -0.00966798 +471 -0.046954 +472 0.0918711 +473 0.0272196 +474 0.0185346 +475 -0.0329217 +476 -0.0536179 +477 0.0803566 +478 -0.0499015 +479 -0.039256 +480 0.0515844 +481 -0.0254139 +482 0.0432167 +483 0.0460773 +484 -0.0744449 +485 0.134048 +486 -0.0109604 +487 -0.000548119 +488 -0.0285135 +489 0.081197 +490 -0.0220405 +491 0.0489457 +492 0.0368539 +493 0.0850633 +494 0.00580247 +495 0.00387908 +496 -0.046728 +497 0.144958 +498 0.08649 +499 -0.0219127 +500 -0.0263977 diff --git a/src/gromacs/gmxana/gmx_analyze.c b/src/gromacs/gmxana/gmx_analyze.c index 02742476b1..c4096954a2 100644 --- a/src/gromacs/gmxana/gmx_analyze.c +++ b/src/gromacs/gmxana/gmx_analyze.c @@ -60,6 +60,7 @@ #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/futil.h" #include "gromacs/utility/smalloc.h" +#include "gromacs/utility/snprintf.h" /* must correspond to char *avbar_opt[] declared in main() */ enum { @@ -394,9 +395,21 @@ static void average(const char *avfile, int avbar_opt, } } -static real anal_ee_inf(double *parm, real T) +/*! \brief Compute final error estimate. + * + * See Eqn A17, Hess, JCP 116 (2002) 209-217 for details. + */ +static real optimal_error_estimate(double sigma, double fitparm[], real tTotal) { - return sqrt(parm[1]*2*parm[0]/T+parm[3]*2*parm[2]/T); + double ss = fitparm[1]*fitparm[0]+(1-fitparm[1])*fitparm[2]; + if ((tTotal <= 0) || (ss <= 0)) + { + fprintf(stderr, "Problem in error estimate: T = %g, ss = %g\n", + tTotal, ss); + return 0; + } + + return sigma*sqrt(2*ss/tTotal); } static void estimate_error(const char *eefile, int nb_min, int resol, int n, @@ -411,7 +424,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, double blav, var; char **leg; real *tbs, *ybs, rtmp, dens, *fitsig, twooe, tau1_est, tau_sig; - double fitparm[4]; + double fitparm[3]; real ee, a, tau1, tau2; if (n < 4) @@ -553,8 +566,8 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, */ fitparm[2] = sqrt(tau1_est*(n-1)*dt); do_lmfit(nbs, ybs, fitsig, 0, tbs, 0, dt*n, oenv, - bDebugMode(), effnERREST, fitparm, 0); - fitparm[3] = 1-fitparm[1]; + bDebugMode(), effnERREST, fitparm, 0, + NULL); } if (bSingleExpFit || fitparm[0] < 0 || fitparm[2] < 0 || fitparm[1] < 0 || (fitparm[1] > 1 && !bAllowNegLTCorr) || fitparm[2] > (n-1)*dt) @@ -573,7 +586,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, fprintf(stdout, "a fitted parameter is negative\n"); } fprintf(stdout, "invalid fit: e.e. %g a %g tau1 %g tau2 %g\n", - sig[s]*anal_ee_inf(fitparm, n*dt), + optimal_error_estimate(sig[s], fitparm, n*dt), fitparm[1], fitparm[0], fitparm[2]); /* Do a fit with tau2 fixed at the total time. * One could also choose any other large value for tau2. @@ -581,10 +594,9 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, fitparm[0] = tau1_est; fitparm[1] = 0.95; fitparm[2] = (n-1)*dt; - fprintf(stderr, "Will fix tau2 at the total time: %g\n", fitparm[2]); + fprintf(stdout, "Will fix tau2 at the total time: %g\n", fitparm[2]); do_lmfit(nbs, ybs, fitsig, 0, tbs, 0, dt*n, oenv, bDebugMode(), - effnERREST, fitparm, 4); - fitparm[3] = 1-fitparm[1]; + effnERREST, fitparm, 4, NULL); } if (bSingleExpFit || fitparm[0] < 0 || fitparm[1] < 0 || (fitparm[1] > 1 && !bAllowNegLTCorr)) @@ -593,7 +605,7 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, { fprintf(stdout, "a fitted parameter is negative\n"); fprintf(stdout, "invalid fit: e.e. %g a %g tau1 %g tau2 %g\n", - sig[s]*anal_ee_inf(fitparm, n*dt), + optimal_error_estimate(sig[s], fitparm, n*dt), fitparm[1], fitparm[0], fitparm[2]); } /* Do a single exponential fit */ @@ -602,11 +614,10 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, fitparm[1] = 1.0; fitparm[2] = 0.0; do_lmfit(nbs, ybs, fitsig, 0, tbs, 0, dt*n, oenv, bDebugMode(), - effnERREST, fitparm, 6); - fitparm[3] = 1-fitparm[1]; + effnERREST, fitparm, 6, NULL); } } - ee = sig[s]*anal_ee_inf(fitparm, n*dt); + ee = optimal_error_estimate(sig[s], fitparm, n*dt); a = fitparm[1]; tau1 = fitparm[0]; tau2 = fitparm[2]; @@ -616,12 +627,14 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, if (output_env_get_xvg_format(oenv) == exvgXMGR) { fprintf(fp, "@ legend string %d \"av %f\"\n", 2*s, av[s]); - fprintf(fp, "@ legend string %d \"ee %6g\"\n", 2*s+1, sig[s]*anal_ee_inf(fitparm, n*dt)); + fprintf(fp, "@ legend string %d \"ee %6g\"\n", 2*s+1, + optimal_error_estimate(sig[s], fitparm, n*dt)); } else if (output_env_get_xvg_format(oenv) == exvgXMGRACE) { fprintf(fp, "@ s%d legend \"av %f\"\n", 2*s, av[s]); - fprintf(fp, "@ s%d legend \"ee %6g\"\n", 2*s+1, sig[s]*anal_ee_inf(fitparm, n*dt)); + fprintf(fp, "@ s%d legend \"ee %6g\"\n", 2*s+1, + optimal_error_estimate(sig[s], fitparm, n*dt)); } for (i = 0; i < nbs; i++) { @@ -672,11 +685,11 @@ static void estimate_error(const char *eefile, int nb_min, int resol, int n, ac_fit[1] = 0.95; ac_fit[2] = 10*acint; do_lmfit(n/nb_min, ac, fitsig, dt, 0, 0, fitlen*dt, oenv, - bDebugMode(), effnEXP3, ac_fit, 0); + bDebugMode(), effnEXPEXP, ac_fit, 0, NULL); ac_fit[3] = 1 - ac_fit[1]; fprintf(stdout, "Set %3d: ac erest %g a %g tau1 %g tau2 %g\n", - s+1, sig[s]*anal_ee_inf(ac_fit, n*dt), + s+1, optimal_error_estimate(sig[s], ac_fit, n*dt), ac_fit[1], ac_fit[0], ac_fit[2]); fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : ""); @@ -789,7 +802,8 @@ static void filter(real flen, int n, int nset, real **val, real dt) static void do_fit(FILE *out, int n, gmx_bool bYdy, int ny, real *x0, real **val, - int npargs, t_pargs *ppa, const output_env_t oenv) + int npargs, t_pargs *ppa, const output_env_t oenv, + const char *fn_fitted) { real *c1 = NULL, *sig = NULL; double *fitparm; @@ -838,7 +852,7 @@ static void do_fit(FILE *out, int n, gmx_bool bYdy, fitparm[0] = 0.5; fitparm[1] = c1[0]; break; - case effnEXP3: + case effnEXPEXP: fitparm[0] = 1.0; fitparm[1] = 0.5*c1[0]; fitparm[2] = 10.0; @@ -877,7 +891,8 @@ static void do_fit(FILE *out, int n, gmx_bool bYdy, fprintf(out, "a%-2d = %12.5e\n", i+1, fitparm[i]); } if (do_lmfit(ny, c1, sig, 0, x0, tbeginfit, tendfit, - oenv, bDebugMode(), efitfn, fitparm, 0)) + oenv, bDebugMode(), efitfn, fitparm, 0, + fn_fitted) > 0) { for (i = 0; (i < nparm); i++) { @@ -1006,6 +1021,50 @@ static void do_geminate(const char *gemFile, int nData, xvgrclose(fp); } +static void print_fitted_function(const char *fitfile, + const char *fn_fitted, + gmx_bool bXYdy, + int nset, + int n, + real *t, + real **val, + int npargs, + t_pargs *ppa, + output_env_t oenv) +{ + FILE *out_fit = gmx_ffopen(fitfile, "w"); + if (bXYdy && nset >= 2) + { + do_fit(out_fit, 0, TRUE, n, t, val, npargs, ppa, oenv, + fn_fitted); + } + else + { + char *buf2 = NULL; + int s, buflen = 0; + if (NULL != fn_fitted) + { + buflen = strlen(fn_fitted)+32; + snew(buf2, buflen); + strncpy(buf2, fn_fitted, buflen); + buf2[strlen(buf2)-4] = '\0'; + } + for (s = 0; s < nset; s++) + { + char *buf = NULL; + if (NULL != fn_fitted) + { + snew(buf, buflen); + snprintf(buf, n, "%s_%d.xvg", buf2, s); + } + do_fit(out_fit, s, FALSE, n, t, val, npargs, ppa, oenv, buf); + sfree(buf); + } + sfree(buf2); + } + gmx_ffclose(out_fit); +} + int gmx_analyze(int argc, char *argv[]) { static const char *desc[] = { @@ -1101,7 +1160,13 @@ int gmx_analyze(int argc, char *argv[]) "Option [TT]-luzar[tt] performs a Luzar & Chandler kinetics analysis", "on output from [gmx-hbond]. The input file can be taken directly", - "from [TT]gmx hbond -ac[tt], and then the same result should be produced." + "from [TT]gmx hbond -ac[tt], and then the same result should be produced.[PAR]", + "Option [TT]-fitfn[tt] performs curve fitting to a number of different", + "curves that make sense in the context of molecular dynamics, mainly", + "exponential curves. More information is in the manual. To check the output", + "of the fitting procedure the option [TT]-fitted[tt] will print both the", + "original data and the fitted function to a new data file. The fitting", + "parameters are stored as comment in the output file." }; static real tb = -1, te = -1, frac = 0.5, filtlen = 0, binwidth = 0.1, aver_start = 0; static gmx_bool bHaveT = TRUE, bDer = FALSE, bSubAv = TRUE, bAverCorr = FALSE, bXYdy = FALSE; @@ -1203,6 +1268,7 @@ int gmx_analyze(int argc, char *argv[]) { efXVG, "-av", "average", ffOPTWR }, { efXVG, "-ee", "errest", ffOPTWR }, { efXVG, "-bal", "ballisitc", ffOPTWR }, + { efXVG, "-fitted", "fitted", ffOPTWR }, /* { efXVG, "-gem", "geminate", ffOPTWR }, */ { efLOG, "-g", "fitlog", ffOPTWR } }; @@ -1283,19 +1349,12 @@ int gmx_analyze(int argc, char *argv[]) if (fitfile != NULL) { - out_fit = gmx_ffopen(fitfile, "w"); - if (bXYdy && nset >= 2) - { - do_fit(out_fit, 0, TRUE, n, t, val, npargs, ppa, oenv); - } - else - { - for (s = 0; s < nset; s++) - { - do_fit(out_fit, s, FALSE, n, t, val, npargs, ppa, oenv); - } - } - gmx_ffclose(out_fit); + print_fitted_function(fitfile, + opt2fn_null("-fitted", NFILE, fnm), + bXYdy, nset, + n, t, val, + npargs, ppa, + oenv); } printf(" std. dev. relative deviation of\n"); diff --git a/src/gromacs/gmxana/gmx_densorder.cpp b/src/gromacs/gmxana/gmx_densorder.cpp index bf4f6c68ea..979f36bc20 100644 --- a/src/gromacs/gmxana/gmx_densorder.cpp +++ b/src/gromacs/gmxana/gmx_densorder.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2010,2011,2012,2013,2014, by the GROMACS development team, led by + * Copyright (c) 2010,2011,2012,2013,2014,2015, 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. @@ -492,9 +492,9 @@ static void interfaces_txy (real ****Densmap, int xslices, int yslices, int zsli /*Fit average density in z over whole trajectory to obtain tentative fit-parameters in fit1 and fit2*/ /*Fit 1st half of box*/ - do_lmfit(zslices, zDensavg, sigma1, binwidth, NULL, startpoint, splitpoint, oenv, FALSE, effnERF, beginfit1, 3); + do_lmfit(zslices, zDensavg, sigma1, binwidth, NULL, startpoint, splitpoint, oenv, FALSE, effnERF, beginfit1, 8, NULL); /*Fit 2nd half of box*/ - do_lmfit(zslices, zDensavg, sigma2, binwidth, NULL, splitpoint, endpoint, oenv, FALSE, effnERF, beginfit2, 3); + do_lmfit(zslices, zDensavg, sigma2, binwidth, NULL, splitpoint, endpoint, oenv, FALSE, effnERF, beginfit2, 8, NULL); /*Initialise the const arrays for storing the average fit parameters*/ avgfit1 = beginfit1; @@ -518,10 +518,10 @@ static void interfaces_txy (real ****Densmap, int xslices, int yslices, int zsli fit2[k] = avgfit2[k]; } /*Now fit and store in structures in row-major order int[n][i][j]*/ - do_lmfit(zslices, Densmap[n][i][j], sigma1, binwidth, NULL, startpoint, splitpoint, oenv, FALSE, effnERF, fit1, 1); + do_lmfit(zslices, Densmap[n][i][j], sigma1, binwidth, NULL, startpoint, splitpoint, oenv, FALSE, effnERF, fit1, 0, NULL); int1[n][j+(yslices*i)]->Z = fit1[2]; int1[n][j+(yslices*i)]->t = fit1[3]; - do_lmfit(zslices, Densmap[n][i][j], sigma2, binwidth, NULL, splitpoint, endpoint, oenv, FALSE, effnERF, fit2, 2); + do_lmfit(zslices, Densmap[n][i][j], sigma2, binwidth, NULL, splitpoint, endpoint, oenv, FALSE, effnERF, fit2, 0, NULL); int2[n][j+(yslices*i)]->Z = fit2[2]; int2[n][j+(yslices*i)]->t = fit2[3]; } diff --git a/src/gromacs/gmxana/gmx_dielectric.c b/src/gromacs/gmxana/gmx_dielectric.c index 7bc7febc79..6e76ed186c 100644 --- a/src/gromacs/gmxana/gmx_dielectric.c +++ b/src/gromacs/gmxana/gmx_dielectric.c @@ -397,7 +397,7 @@ int gmx_dielectric(int argc, char *argv[]) integral = print_and_integrate(NULL, calc_nbegin(nx, y[0], tbegin), dt, y[1], NULL, 1); integral += do_lmfit(nx, y[1], y[2], dt, y[0], tbegin, tend, - oenv, TRUE, eFitFn, fitparms, fix); + oenv, TRUE, eFitFn, fitparms, fix, NULL); for (i = 0; i < nx; i++) { y[3][i] = fit_function(eFitFn, fitparms, y[0][i]); diff --git a/src/gromacs/gmxana/gmx_hbond.c b/src/gromacs/gmxana/gmx_hbond.c index f9e8fb89c6..0d364128b1 100644 --- a/src/gromacs/gmxana/gmx_hbond.c +++ b/src/gromacs/gmxana/gmx_hbond.c @@ -2402,7 +2402,7 @@ static void smooth_tail(int n, real t[], real c[], real sigma_c[], real start, } fitparm[1] = 0.95; do_lmfit(n, c, sigma_c, 0, t, start, t[n-1], oenv, bDebugMode(), - effnEXP2, fitparm, 0); + effnEXP2, fitparm, 0, NULL); } void analyse_corr(int n, real t[], real ct[], real nt[], real kt[], diff --git a/src/gromacs/gmxana/gmx_tcaf.c b/src/gromacs/gmxana/gmx_tcaf.c index d2a6b72d74..eef238878a 100644 --- a/src/gromacs/gmxana/gmx_tcaf.c +++ b/src/gromacs/gmxana/gmx_tcaf.c @@ -208,7 +208,7 @@ static void process_tcaf(int nframes, real dt, int nkc, real **tc, rvec *kfac, fitparms[0] = 1; fitparms[1] = 1; do_lmfit(ncorr, tcaf[k], sig, dt, 0, 0, ncorr*dt, - oenv, bDebugMode(), effnVAC, fitparms, 0); + oenv, bDebugMode(), effnVAC, fitparms, 0, NULL); eta = 1000*fitparms[1]*rho/ (4*fitparms[0]*PICO*norm2(kfac[k])/(NANO*NANO)); fprintf(stdout, "k %6.3f tau %6.3f eta %8.5f 10^-3 kg/(m s)\n", @@ -233,7 +233,7 @@ static void process_tcaf(int nframes, real dt, int nkc, real **tc, rvec *kfac, fitparms[0] = 1; fitparms[1] = 1; do_lmfit(ncorr, tcafc[k], sig, dt, 0, 0, ncorr*dt, - oenv, bDebugMode(), effnVAC, fitparms, 0); + oenv, bDebugMode(), effnVAC, fitparms, 0, NULL); eta = 1000*fitparms[1]*rho/ (4*fitparms[0]*PICO*norm2(kfac[kset_c[k]])/(NANO*NANO)); fprintf(stdout, -- 2.22.0