src/gromacs/correlationfunctions/integrate.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2017,2018, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \brief
  38  * Implement routines for integrating a data set
  39  *
  40  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "integrate.h"
  45
  46 #include <cmath>
  47 #include <cstdio>
  48
  49 #include "gromacs/math/functions.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/utility/fatalerror.h"
  52
  53 /*! \brief Integrate a function and printe the integral value. */
  54 real print_and_integrate(FILE *fp, int n, real dt, const real c[],
  55                          const real *fit, int nskip)
  56 {
  57     real c0, sum;
  58     int  j;
  59
  60     /* Use trapezoidal rule for calculating integral */
  61     sum = 0.0;
  62     for (j = 0; (j < n); j++)
  63     {
  64         c0 = c[j];
  65         if (fp && (nskip == 0 || j % nskip == 0))
  66         {
  67             fprintf(fp, "%10.3f  %10.5f\n", j*dt, c0);
  68         }
  69         if (j > 0)
  70         {
  71             sum += dt*(c0+c[j-1]);
  72         }
  73     }
  74     if (fp)
  75     {
  76         fprintf(fp, "&\n");
  77         if (fit)
  78         {
  79             for (j = 0; (j < n); j++)
  80             {
  81                 if (nskip == 0 || j % nskip == 0)
  82                 {
  83                     fprintf(fp, "%10.3f  %10.5f\n", j*dt, fit[j]);
  84                 }
  85             }
  86             fprintf(fp, "&\n");
  87         }
  88     }
  89     return sum*0.5;
  90 }
  91
  92 /*! \brief Compute and return the integral of a function. */
  93 real evaluate_integral(int n, const real x[], const real y[],
  94                        const real dy[], real aver_start,
  95                        real *stddev)
  96 {
  97     double sum, sum_var, w;
  98     double sum_tail = 0, sum2_tail = 0;
  99     int    j, nsum_tail = 0;
 100
 101     /* Use trapezoidal rule for calculating integral */
 102     if (n <= 0)
 103     {
 104         gmx_fatal(FARGS, "Evaluating integral: n = %d (file %s, line %d)",
 105                   n, __FILE__, __LINE__);
 106     }
 107
 108     sum     = 0;
 109     sum_var = 0;
 110     for (j = 0; (j < n); j++)
 111     {
 112         w = 0;
 113         if (j > 0)
 114         {
 115             w += 0.5*(x[j] - x[j-1]);
 116         }
 117         if (j < n-1)
 118         {
 119             w += 0.5*(x[j+1] - x[j]);
 120         }
 121         sum += w*y[j];
 122         if (dy)
 123         {
 124             /* Assume all errors are uncorrelated */
 125             sum_var += gmx::square(w*dy[j]);
 126         }
 127
 128         if ((aver_start > 0) && (x[j] >= aver_start))
 129         {
 130             sum_tail  += sum;
 131             sum2_tail += std::sqrt(sum_var);
 132             nsum_tail += 1;
 133         }
 134     }
 135
 136     if (nsum_tail > 0)
 137     {
 138         sum = sum_tail/nsum_tail;
 139         /* This is a worst case estimate, assuming all stddev's are correlated. */
 140         *stddev = sum2_tail/nsum_tail;
 141     }
 142     else
 143     {
 144         *stddev = std::sqrt(sum_var);
 145     }
 146
 147     return sum;
 148 }