Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_bar.c
index 45ae72b2242623fccd0b39512afc78df028d8ac4..2a2b8fd8cb212cfaa1f36a9d373cbf992afd0fac 100644 (file)
@@ -1,64 +1,60 @@
-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * This file is part of the GROMACS molecular simulation package.
  *
+ * Copyright (c) 2010,2011,2012,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.
  *
- *                This source code is part of
- *
- *                 G   R   O   M   A   C   S
- *
- *          GROningen MAchine for Chemical Simulations
- *
- *                        VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * 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.
  *
- * If you want to redistribute modifications, 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 www.gromacs.org.
+ * 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.
  *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
+ * 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.
  *
- * For more info, check our website at http://www.gromacs.org
+ * 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.
  *
- * And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
 
-#include <math.h>
-#include <string.h>
 #include <ctype.h>
-#include <math.h>
 #include <float.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
 
-#include "sysstuff.h"
-#include "typedefs.h"
-#include "smalloc.h"
-#include "futil.h"
-#include "statutil.h"
-#include "copyrite.h"
-#include "macros.h"
-#include "enxio.h"
-#include "physics.h"
-#include "gmx_fatal.h"
-#include "xvgr.h"
-#include "gmx_ana.h"
-#include "maths.h"
-#include "string2.h"
-#include "names.h"
-#include "mdebin.h"
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/fileio/enxio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/mdebin.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/viewit.h"
+#include "gromacs/math/units.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
 
 
 /* Structure for the names of lambda vector components */
@@ -108,12 +104,12 @@ typedef struct hist_t
     unsigned int   *bin[2];                 /* the (forward + reverse) histogram values */
     double          dx[2];                  /* the histogram spacing. The reverse
                                                dx is the negative of the forward dx.*/
-    gmx_large_int_t x0[2];                  /* the (forward + reverse) histogram start
-                                               point(s) as int */
+    gmx_int64_t     x0[2];                  /* the (forward + reverse) histogram start
+                                                   point(s) as int */
 
     int             nbin[2];                /* the (forward+reverse) number of bins */
-    gmx_large_int_t sum;                    /* the total number of counts. Must be
-                                               the same for forward + reverse.  */
+    gmx_int64_t     sum;                    /* the total number of counts. Must be
+                                                   the same for forward + reverse.  */
     int             nhist;                  /* number of hist datas (forward or reverse) */
 
     double          start_time, delta_time; /* start time, end time of histogram */
@@ -142,7 +138,7 @@ typedef struct samples_t
     size_t          ndu_alloc, nt_alloc; /* pre-allocated sizes */
     hist_t         *hist_alloc;          /* allocated hist */
 
-    gmx_large_int_t ntot;                /* total number of samples */
+    gmx_int64_t     ntot;                /* total number of samples */
     const char     *filename;            /* the file name this sample comes from */
 } samples_t;
 
@@ -172,8 +168,8 @@ typedef struct sample_coll_t
     sample_range_t *r;                 /* the sample ranges */
     int             nsamples_alloc;    /* number of allocated samples */
 
-    gmx_large_int_t ntot;              /* total number of samples in the ranges of
-                                          this collection */
+    gmx_int64_t     ntot;              /* total number of samples in the ranges of
+                                              this collection */
 
     struct sample_coll_t *next, *prev; /* next and previous in the list */
 } sample_coll_t;
@@ -236,7 +232,7 @@ static void lambda_components_add(lambda_components_t *lc,
     while (lc->N + 1 > lc->Nalloc)
     {
         lc->Nalloc = (lc->Nalloc == 0) ? 2 : 2*lc->Nalloc;
-        srealloc( lc->names, lc->Nalloc );
+        srenew(lc->names, lc->Nalloc);
     }
     snew(lc->names[lc->N], name_length+1);
     strncpy(lc->names[lc->N], name, name_length);
@@ -358,7 +354,7 @@ static void lambda_vec_print(const lambda_vec_t *lv, char *str, gmx_bool named)
         str += sprintf(str, "dH/dl");
         if (strlen(lv->lc->names[lv->dhdl]) > 0)
         {
-            str += sprintf(str, " (%s)", lv->lc->names[lv->dhdl]);
+            sprintf(str, " (%s)", lv->lc->names[lv->dhdl]);
         }
     }
 }
@@ -1134,6 +1130,31 @@ void sim_data_histogram(sim_data_t *sd, const char *filename,
     xvgrclose(fp);
 }
 
+static int
+snprint_lambda_vec(char *str, int sz, const char *label, lambda_vec_t *lambda)
+{
+    int n = 0;
+
+    n += snprintf(str+n, sz-n, "lambda vector [%s]: ", label);
+    if (lambda->index >= 0)
+    {
+        n += snprintf(str+n, sz-n, " init-lambda-state=%d", lambda->index);
+    }
+    if (lambda->dhdl >= 0)
+    {
+        n += snprintf(str+n, sz-n, " dhdl index=%d", lambda->dhdl);
+    }
+    else
+    {
+        int i;
+        for (i = 0; i < lambda->lc->N; i++)
+        {
+            n += snprintf(str+n, sz-n, " (%s) l=%g", lambda->lc->names[i], lambda->val[i]);
+        }
+    }
+    return n;
+}
+
 /* create a collection (array) of barres_t object given a ordered linked list
    of barlamda_t sample collections */
 static barres_t *barres_list_create(sim_data_t *sd, int *nres,
@@ -1189,17 +1210,27 @@ static barres_t *barres_list_create(sim_data_t *sd, int *nres,
         }
         else if (!scprev && !sc)
         {
-            gmx_fatal(FARGS, "There is no path from lambda=%g -> %g that is covered by foreign lambdas:\ncannot proceed with BAR.\nUse thermodynamic integration of dH/dl by calculating the averages of dH/dl\nwith g_analyze and integrating them.\nAlternatively, use the -extp option if (and only if) the Hamiltonian\ndepends linearly on lambda, which is NOT normally the case.\n", bl->prev->lambda, bl->lambda);
+            char descX[STRLEN], descY[STRLEN];
+            snprint_lambda_vec(descX, STRLEN, "X", bl->prev->lambda);
+            snprint_lambda_vec(descY, STRLEN, "Y", bl->lambda);
+
+            gmx_fatal(FARGS, "There is no path between the states X & Y below that is covered by foreign lambdas:\ncannot proceed with BAR.\nUse thermodynamic integration of dH/dl by calculating the averages of dH/dl\nwith g_analyze and integrating them.\nAlternatively, use the -extp option if (and only if) the Hamiltonian\ndepends linearly on lambda, which is NOT normally the case.\n\n%s\n%s\n", descX, descY);
         }
 
         /* normal delta H */
         if (!scprev)
         {
-            gmx_fatal(FARGS, "Could not find a set for foreign lambda = %g\nin the files for lambda = %g", bl->lambda, bl->prev->lambda);
+            char descX[STRLEN], descY[STRLEN];
+            snprint_lambda_vec(descX, STRLEN, "X", bl->lambda);
+            snprint_lambda_vec(descY, STRLEN, "Y", bl->prev->lambda);
+            gmx_fatal(FARGS, "Could not find a set for foreign lambda (state X below)\nin the files for main lambda (state Y below)\n\n%s\n%s\n", descX, descY);
         }
         if (!sc)
         {
-            gmx_fatal(FARGS, "Could not find a set for foreign lambda = %g\nin the files for lambda = %g", bl->prev->lambda, bl->lambda);
+            char descX[STRLEN], descY[STRLEN];
+            snprint_lambda_vec(descX, STRLEN, "X", bl->prev->lambda);
+            snprint_lambda_vec(descY, STRLEN, "Y", bl->lambda);
+            gmx_fatal(FARGS, "Could not find a set for foreign lambda (state X below)\nin the files for main lambda (state Y below)\n\n%s\n%s\n", descX, descY);
         }
         br->a = scprev;
         br->b = sc;
@@ -1421,9 +1452,9 @@ static gmx_bool sample_coll_create_subsample(sample_coll_t  *sc,
     int             j;
     int             hist_start, hist_end;
 
-    gmx_large_int_t ntot_start;
-    gmx_large_int_t ntot_end;
-    gmx_large_int_t ntot_so_far;
+    gmx_int64_t     ntot_start;
+    gmx_int64_t     ntot_end;
+    gmx_int64_t     ntot_so_far;
 
     *sc = *sc_orig; /* just copy all fields */
 
@@ -1440,13 +1471,13 @@ static gmx_bool sample_coll_create_subsample(sample_coll_t  *sc,
 
     /* now fix start and end fields */
     /* the casts avoid possible overflows */
-    ntot_start  = (gmx_large_int_t)(sc_orig->ntot*(double)i/(double)ni);
-    ntot_end    = (gmx_large_int_t)(sc_orig->ntot*(double)(i+1)/(double)ni);
+    ntot_start  = (gmx_int64_t)(sc_orig->ntot*(double)i/(double)ni);
+    ntot_end    = (gmx_int64_t)(sc_orig->ntot*(double)(i+1)/(double)ni);
     ntot_so_far = 0;
     for (j = 0; j < sc->nsamples; j++)
     {
-        gmx_large_int_t ntot_add;
-        gmx_large_int_t new_start, new_end;
+        gmx_int64_t ntot_add;
+        gmx_int64_t new_start, new_end;
 
         if (sc->r[j].use)
         {
@@ -1733,9 +1764,6 @@ static double calc_bar_lowlevel(sample_coll_t *ca, sample_coll_t *cb,
     {
         DG1 = 0.5*(DG0 + DG2);
 
-        /*printf("Wfac1=%g, Wfac2=%g, beta=%g, DG1=%g\n",Wfac1,Wfac2,beta,
-           DG1);*/
-
         /* calculate the BAR averages */
         dDG1 = 0.;
 
@@ -2437,7 +2465,6 @@ static gmx_bool read_lambda_vector(const char   *str,
  */
 static gmx_bool legend2lambda(const char   *fn,
                               const char   *legend,
-                              xvg_t        *ba,
                               lambda_vec_t *lam)
 {
     double        lambda = 0;
@@ -2683,7 +2710,7 @@ static gmx_bool subtitle2lambda(const char *subtitle, xvg_t *ba, const char *fn,
     return bFound;
 }
 
-static void filename2lambda(const char *fn, xvg_t *ba)
+static void filename2lambda(const char *fn)
 {
     double        lambda;
     const char   *ptr, *digitptr;
@@ -2779,7 +2806,7 @@ static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
             {
                 if (ba->temp <= 0)
                 {
-                    gmx_fatal(FARGS, "Found temperature of %g in file '%s'",
+                    gmx_fatal(FARGS, "Found temperature of %f in file '%s'",
                               ba->temp, fn);
                 }
             }
@@ -2789,7 +2816,7 @@ static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
     {
         if (*temp <= 0)
         {
-            gmx_fatal(FARGS, "Did not find a temperature in the subtitle in file '%s', use the -temp option of [TT]g_bar[tt]", fn);
+            gmx_fatal(FARGS, "Did not find a temperature in the subtitle in file '%s', use the -temp option of [TT]gmx bar[tt]", fn);
         }
         ba->temp = *temp;
     }
@@ -2811,7 +2838,7 @@ static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
             if (!native_lambda_read)
             {
                 /* Deduce lambda from the file name */
-                filename2lambda(fn, ba);
+                filename2lambda(fn);
                 native_lambda_read = TRUE;
             }
             ba->lambda[0] = ba->native_lambda;
@@ -2829,7 +2856,7 @@ static void read_bar_xvg_lowlevel(const char *fn, real *temp, xvg_t *ba,
             /* Read lambda from the legend */
             lambda_vec_init( &(ba->lambda[i]), lc );
             lambda_vec_copy( &(ba->lambda[i]), &(ba->native_lambda));
-            use = legend2lambda(fn, legend[i], ba, &(ba->lambda[i]));
+            use = legend2lambda(fn, legend[i], &(ba->lambda[i]));
             if (use)
             {
                 lambda_vec_print(&(ba->lambda[i]), buf, FALSE);
@@ -2940,7 +2967,7 @@ static void read_edr_rawdh_block(samples_t **smp, int *ndu, t_enxblock *blk,
          (blk->sub[1].nr < 1) )
     {
         gmx_fatal(FARGS,
-                  "Unexpected/corrupted block data in file %s around time %g.",
+                  "Unexpected/corrupted block data in file %s around time %f.",
                   filename, start_time);
     }
 
@@ -2987,7 +3014,7 @@ static void read_edr_rawdh_block(samples_t **smp, int *ndu, t_enxblock *blk,
         lambda_vec_print(foreign_lambda, buf, FALSE);
         lambda_vec_print(s->foreign_lambda, buf2, FALSE);
         fprintf(stderr, "Got foreign lambda=%s, expected: %s\n", buf, buf2);
-        gmx_fatal(FARGS, "Corrupted data in file %s around t=%g.",
+        gmx_fatal(FARGS, "Corrupted data in file %s around t=%f.",
                   filename, start_time);
     }
 
@@ -3038,12 +3065,12 @@ static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
     /* check the block types etc. */
     if ( (blk->nsub < 2) ||
          (blk->sub[0].type != xdr_datatype_double) ||
-         (blk->sub[1].type != xdr_datatype_large_int) ||
+         (blk->sub[1].type != xdr_datatype_int64) ||
          (blk->sub[0].nr < 2)  ||
          (blk->sub[1].nr < 2) )
     {
         gmx_fatal(FARGS,
-                  "Unexpected/corrupted block data in file %s around time %g",
+                  "Unexpected/corrupted block data in file %s around time %f",
                   filename, start_time);
     }
 
@@ -3055,7 +3082,7 @@ static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
     if (nhist > 2)
     {
         gmx_fatal(FARGS,
-                  "Unexpected/corrupted block data in file %s around time %g",
+                  "Unexpected/corrupted block data in file %s around time %f",
                   filename, start_time);
     }
 
@@ -3136,7 +3163,7 @@ static samples_t *read_edr_hist_block(int *nsamples, t_enxblock *blk,
     for (i = 0; i < nhist; i++)
     {
         int             nbin;
-        gmx_large_int_t sum = 0;
+        gmx_int64_t     sum = 0;
 
         for (j = 0; j < s->hist->nbin[i]; j++)
         {
@@ -3272,7 +3299,7 @@ static void read_barsim_edr(char *fn, real *temp, sim_data_t *sd)
                     if (fr->block[i].nsub < 2)
                     {
                         gmx_fatal(FARGS,
-                                  "No lambda vector, but start_lambda=%g\n",
+                                  "No lambda vector, but start_lambda=%f\n",
                                   old_start_lambda);
                     }
                     n_lambda_vec = fr->block[i].sub[1].ival[1];
@@ -3320,7 +3347,7 @@ static void read_barsim_edr(char *fn, real *temp, sim_data_t *sd)
             /* check the native lambda */
             if (!lambda_vec_same(&start_lambda, native_lambda) )
             {
-                gmx_fatal(FARGS, "Native lambda not constant in file %s: started at %g, and becomes %g at time %g",
+                gmx_fatal(FARGS, "Native lambda not constant in file %s: started at %f, and becomes %f at time %f",
                           fn, native_lambda, start_lambda, start_time);
             }
             /* check the number of samples against the previous number */
@@ -3464,7 +3491,7 @@ static void read_barsim_edr(char *fn, real *temp, sim_data_t *sd)
 int gmx_bar(int argc, char *argv[])
 {
     static const char *desc[] = {
-        "[TT]g_bar[tt] calculates free energy difference estimates through ",
+        "[THISMODULE] calculates free energy difference estimates through ",
         "Bennett's acceptance ratio method (BAR). It also automatically",
         "adds series of individual free energies obtained with BAR into",
         "a combined free energy estimate.[PAR]",
@@ -3526,7 +3553,7 @@ int gmx_bar(int argc, char *argv[])
         "over 5 blocks. A range of block numbers for error estimation can ",
         "be provided with the options [TT]-nbmin[tt] and [TT]-nbmax[tt].[PAR]",
 
-        "[TT]g_bar[tt] tries to aggregate samples with the same 'native' and ",
+        "[THISMODULE] tries to aggregate samples with the same 'native' and ",
         "'foreign' [GRK]lambda[grk] values, but always assumes independent ",
         "samples. [BB]Note[bb] that when aggregating energy ",
         "differences/derivatives with different sampling intervals, this is ",
@@ -3544,7 +3571,7 @@ int gmx_bar(int argc, char *argv[])
         "[TT]*[tt]  lam_B: the [GRK]lambda[grk] values for point B.[BR]",
         "[TT]*[tt]     DG: the free energy estimate.[BR]",
         "[TT]*[tt]    s_A: an estimate of the relative entropy of B in A.[BR]",
-        "[TT]*[tt]    s_A: an estimate of the relative entropy of A in B.[BR]",
+        "[TT]*[tt]    s_B: an estimate of the relative entropy of A in B.[BR]",
         "[TT]*[tt]  stdev: an estimate expected per-sample standard deviation.[PAR]",
 
         "The relative entropy of both states in each other's ensemble can be ",
@@ -3618,9 +3645,12 @@ int gmx_bar(int argc, char *argv[])
     double       sum_histrange_err = 0.; /* histogram range error */
     double       stat_err          = 0.; /* statistical error */
 
-    parse_common_args(&argc, argv,
-                      PCA_CAN_VIEW,
-                      NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
+    if (!parse_common_args(&argc, argv,
+                           PCA_CAN_VIEW,
+                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
+    {
+        return 0;
+    }
 
     if (opt2bSet("-f", NFILE, fnm))
     {
@@ -3941,13 +3971,11 @@ int gmx_bar(int argc, char *argv[])
     {
         lambda_vec_print_short(results[nresults-1].b->native_lambda, buf);
         fprintf(fpi, xvg2format, buf, dg_tot);
-        ffclose(fpi);
+        gmx_ffclose(fpi);
     }
 
     do_view(oenv, opt2fn_null("-o", NFILE, fnm), "-xydy");
     do_view(oenv, opt2fn_null("-oi", NFILE, fnm), "-xydy");
 
-    thanx(stderr);
-
     return 0;
 }