src/gromacs/gmxana/gmx_energy.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, 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 #include "gmxpre.h"
  38
  39 #include <cmath>
  40 #include <cstdlib>
  41 #include <cstring>
  42
  43 #include <algorithm>
  44
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/correlationfunctions/autocorr.h"
  47 #include "gromacs/fileio/enxio.h"
  48 #include "gromacs/fileio/gmxfio.h"
  49 #include "gromacs/fileio/tpxio.h"
  50 #include "gromacs/fileio/trxio.h"
  51 #include "gromacs/fileio/xvgr.h"
  52 #include "gromacs/gmxana/gmx_ana.h"
  53 #include "gromacs/gmxana/gstat.h"
  54 #include "gromacs/legacyheaders/copyrite.h"
  55 #include "gromacs/legacyheaders/macros.h"
  56 #include "gromacs/legacyheaders/mdebin.h"
  57 #include "gromacs/legacyheaders/names.h"
  58 #include "gromacs/legacyheaders/typedefs.h"
  59 #include "gromacs/legacyheaders/viewit.h"
  60 #include "gromacs/math/units.h"
  61 #include "gromacs/math/vec.h"
  62 #include "gromacs/topology/mtop_util.h"
  63 #include "gromacs/utility/cstringutil.h"
  64 #include "gromacs/utility/fatalerror.h"
  65 #include "gromacs/utility/gmxassert.h"
  66 #include "gromacs/utility/smalloc.h"
  67
  68 static real       minthird = -1.0/3.0, minsixth = -1.0/6.0;
  69
  70 typedef struct {
  71     real sum;
  72     real sum2;
  73 } exactsum_t;
  74
  75 typedef struct {
  76     real       *ener;
  77     exactsum_t *es;
  78     gmx_bool    bExactStat;
  79     double      av;
  80     double      rmsd;
  81     double      ee;
  82     double      slope;
  83 } enerdat_t;
  84
  85 typedef struct {
  86     gmx_int64_t      nsteps;
  87     gmx_int64_t      npoints;
  88     int              nframes;
  89     int             *step;
  90     int             *steps;
  91     int             *points;
  92     enerdat_t       *s;
  93     gmx_bool         bHaveSums;
  94 } enerdata_t;
  95
  96 static double mypow(double x, double y)
  97 {
  98     if (x > 0)
  99     {
 100         return std::pow(x, y);
 101     }
 102     else
 103     {
 104         return 0.0;
 105     }
 106 }
 107
 108 static int *select_it(int nre, char *nm[], int *nset)
 109 {
 110     gmx_bool *bE;
 111     int       n, k, j, i;
 112     int      *set;
 113     gmx_bool  bVerbose = TRUE;
 114
 115     if ((getenv("GMX_ENER_VERBOSE")) != NULL)
 116     {
 117         bVerbose = FALSE;
 118     }
 119
 120     fprintf(stderr, "Select the terms you want from the following list\n");
 121     fprintf(stderr, "End your selection with 0\n");
 122
 123     if (bVerbose)
 124     {
 125         for (k = 0; (k < nre); )
 126         {
 127             for (j = 0; (j < 4) && (k < nre); j++, k++)
 128             {
 129                 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
 130             }
 131             fprintf(stderr, "\n");
 132         }
 133     }
 134
 135     snew(bE, nre);
 136     do
 137     {
 138         if (1 != scanf("%d", &n))
 139         {
 140             gmx_fatal(FARGS, "Error reading user input");
 141         }
 142         if ((n > 0) && (n <= nre))
 143         {
 144             bE[n-1] = TRUE;
 145         }
 146     }
 147     while (n != 0);
 148
 149     snew(set, nre);
 150     for (i = (*nset) = 0; (i < nre); i++)
 151     {
 152         if (bE[i])
 153         {
 154             set[(*nset)++] = i;
 155         }
 156     }
 157
 158     sfree(bE);
 159
 160     return set;
 161 }
 162
 163 static void chomp(char *buf)
 164 {
 165     int len = std::strlen(buf);
 166
 167     while ((len > 0) && (buf[len-1] == '\n'))
 168     {
 169         buf[len-1] = '\0';
 170         len--;
 171     }
 172 }
 173
 174 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
 175 {
 176     gmx_bool   *bE;
 177     int         k, kk, j, i, nmatch, nind, nss;
 178     int        *set;
 179     gmx_bool    bEOF, bVerbose = TRUE, bLong = FALSE;
 180     char       *ptr, buf[STRLEN];
 181     const char *fm4   = "%3d  %-14s";
 182     const char *fm2   = "%3d  %-34s";
 183     char      **newnm = NULL;
 184
 185     if ((getenv("GMX_ENER_VERBOSE")) != NULL)
 186     {
 187         bVerbose = FALSE;
 188     }
 189
 190     fprintf(stderr, "\n");
 191     fprintf(stderr, "Select the terms you want from the following list by\n");
 192     fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
 193     fprintf(stderr, "End your selection with an empty line or a zero.\n");
 194     fprintf(stderr, "-------------------------------------------------------------------\n");
 195
 196     snew(newnm, nre);
 197     j = 0;
 198     for (k = 0; k < nre; k++)
 199     {
 200         newnm[k] = gmx_strdup(nm[k].name);
 201         /* Insert dashes in all the names */
 202         while ((ptr = std::strchr(newnm[k], ' ')) != NULL)
 203         {
 204             *ptr = '-';
 205         }
 206         if (bVerbose)
 207         {
 208             if (j == 0)
 209             {
 210                 if (k > 0)
 211                 {
 212                     fprintf(stderr, "\n");
 213                 }
 214                 bLong = FALSE;
 215                 for (kk = k; kk < k+4; kk++)
 216                 {
 217                     if (kk < nre && std::strlen(nm[kk].name) > 14)
 218                     {
 219                         bLong = TRUE;
 220                     }
 221                 }
 222             }
 223             else
 224             {
 225                 fprintf(stderr, " ");
 226             }
 227             if (!bLong)
 228             {
 229                 fprintf(stderr, fm4, k+1, newnm[k]);
 230                 j++;
 231                 if (j == 4)
 232                 {
 233                     j = 0;
 234                 }
 235             }
 236             else
 237             {
 238                 fprintf(stderr, fm2, k+1, newnm[k]);
 239                 j++;
 240                 if (j == 2)
 241                 {
 242                     j = 0;
 243                 }
 244             }
 245         }
 246     }
 247     if (bVerbose)
 248     {
 249         fprintf(stderr, "\n\n");
 250     }
 251
 252     snew(bE, nre);
 253
 254     bEOF = FALSE;
 255     while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
 256     {
 257         /* Remove newlines */
 258         chomp(buf);
 259
 260         /* Remove spaces */
 261         trim(buf);
 262
 263         /* Empty line means end of input */
 264         bEOF = (std::strlen(buf) == 0);
 265         if (!bEOF)
 266         {
 267             ptr = buf;
 268             do
 269             {
 270                 if (!bEOF)
 271                 {
 272                     /* First try to read an integer */
 273                     nss   = sscanf(ptr, "%d", &nind);
 274                     if (nss == 1)
 275                     {
 276                         /* Zero means end of input */
 277                         if (nind == 0)
 278                         {
 279                             bEOF = TRUE;
 280                         }
 281                         else if ((1 <= nind) && (nind <= nre))
 282                         {
 283                             bE[nind-1] = TRUE;
 284                         }
 285                         else
 286                         {
 287                             fprintf(stderr, "number %d is out of range\n", nind);
 288                         }
 289                     }
 290                     else
 291                     {
 292                         /* Now try to read a string */
 293                         nmatch = 0;
 294                         for (nind = 0; nind < nre; nind++)
 295                         {
 296                             if (gmx_strcasecmp(newnm[nind], ptr) == 0)
 297                             {
 298                                 bE[nind] = TRUE;
 299                                 nmatch++;
 300                             }
 301                         }
 302                         if (nmatch == 0)
 303                         {
 304                             i      = std::strlen(ptr);
 305                             nmatch = 0;
 306                             for (nind = 0; nind < nre; nind++)
 307                             {
 308                                 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
 309                                 {
 310                                     bE[nind] = TRUE;
 311                                     nmatch++;
 312                                 }
 313                             }
 314                             if (nmatch == 0)
 315                             {
 316                                 fprintf(stderr, "String '%s' does not match anything\n", ptr);
 317                             }
 318                         }
 319                     }
 320                 }
 321                 /* Look for the first space, and remove spaces from there */
 322                 if ((ptr = std::strchr(ptr, ' ')) != NULL)
 323                 {
 324                     trim(ptr);
 325                 }
 326             }
 327             while (!bEOF && (ptr && (std::strlen(ptr) > 0)));
 328         }
 329     }
 330
 331     snew(set, nre);
 332     for (i = (*nset) = 0; (i < nre); i++)
 333     {
 334         if (bE[i])
 335         {
 336             set[(*nset)++] = i;
 337         }
 338     }
 339
 340     sfree(bE);
 341
 342     if (*nset == 0)
 343     {
 344         gmx_fatal(FARGS, "No energy terms selected");
 345     }
 346
 347     for (i = 0; (i < nre); i++)
 348     {
 349         sfree(newnm[i]);
 350     }
 351     sfree(newnm);
 352
 353     return set;
 354 }
 355
 356 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
 357 {
 358     gmx_mtop_t  mtop;
 359     int         natoms;
 360     matrix      box;
 361
 362     /* all we need is the ir to be able to write the label */
 363     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
 364 }
 365
 366 static void get_orires_parms(const char *topnm,
 367                              int *nor, int *nex, int **label, real **obs)
 368 {
 369     gmx_mtop_t      mtop;
 370     gmx_localtop_t *top;
 371     t_inputrec      ir;
 372     t_iparams      *ip;
 373     int             natoms, i;
 374     t_iatom        *iatom;
 375     int             nb;
 376     matrix          box;
 377
 378     read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
 379     top = gmx_mtop_generate_local_top(&mtop, &ir);
 380
 381     ip       = top->idef.iparams;
 382     iatom    = top->idef.il[F_ORIRES].iatoms;
 383
 384     /* Count how many distance restraint there are... */
 385     nb = top->idef.il[F_ORIRES].nr;
 386     if (nb == 0)
 387     {
 388         gmx_fatal(FARGS, "No orientation restraints in topology!\n");
 389     }
 390
 391     *nor = nb/3;
 392     *nex = 0;
 393     snew(*label, *nor);
 394     snew(*obs, *nor);
 395     for (i = 0; i < nb; i += 3)
 396     {
 397         (*label)[i/3] = ip[iatom[i]].orires.label;
 398         (*obs)[i/3]   = ip[iatom[i]].orires.obs;
 399         if (ip[iatom[i]].orires.ex >= *nex)
 400         {
 401             *nex = ip[iatom[i]].orires.ex+1;
 402         }
 403     }
 404     fprintf(stderr, "Found %d orientation restraints with %d experiments",
 405             *nor, *nex);
 406 }
 407
 408 static int get_bounds(const char *topnm,
 409                       real **bounds, int **index, int **dr_pair, int *npairs,
 410                       gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
 411 {
 412     gmx_localtop_t *top;
 413     t_functype     *functype;
 414     t_iparams      *ip;
 415     int             natoms, i, j, k, type, ftype, natom;
 416     t_ilist        *disres;
 417     t_iatom        *iatom;
 418     real           *b;
 419     int            *ind, *pair;
 420     int             nb, label1;
 421     matrix          box;
 422
 423     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
 424     snew(*ltop, 1);
 425     top   = gmx_mtop_generate_local_top(mtop, ir);
 426     *ltop = top;
 427
 428     functype = top->idef.functype;
 429     ip       = top->idef.iparams;
 430
 431     /* Count how many distance restraint there are... */
 432     nb = top->idef.il[F_DISRES].nr;
 433     if (nb == 0)
 434     {
 435         gmx_fatal(FARGS, "No distance restraints in topology!\n");
 436     }
 437
 438     /* Allocate memory */
 439     snew(b, nb);
 440     snew(ind, nb);
 441     snew(pair, nb+1);
 442
 443     /* Fill the bound array */
 444     nb = 0;
 445     for (i = 0; (i < top->idef.ntypes); i++)
 446     {
 447         ftype = functype[i];
 448         if (ftype == F_DISRES)
 449         {
 450
 451             label1  = ip[i].disres.label;
 452             b[nb]   = ip[i].disres.up1;
 453             ind[nb] = label1;
 454             nb++;
 455         }
 456     }
 457     *bounds = b;
 458
 459     /* Fill the index array */
 460     label1  = -1;
 461     disres  = &(top->idef.il[F_DISRES]);
 462     iatom   = disres->iatoms;
 463     for (i = j = k = 0; (i < disres->nr); )
 464     {
 465         type  = iatom[i];
 466         ftype = top->idef.functype[type];
 467         natom = interaction_function[ftype].nratoms+1;
 468         if (label1 != top->idef.iparams[type].disres.label)
 469         {
 470             pair[j] = k;
 471             label1  = top->idef.iparams[type].disres.label;
 472             j++;
 473         }
 474         k++;
 475         i += natom;
 476     }
 477     pair[j]  = k;
 478     *npairs  = k;
 479     if (j != nb)
 480     {
 481         gmx_incons("get_bounds for distance restraints");
 482     }
 483
 484     *index   = ind;
 485     *dr_pair = pair;
 486
 487     return nb;
 488 }
 489
 490 static void calc_violations(real rt[], real rav3[], int nb, int index[],
 491                             real bounds[], real *viol, double *st, double *sa)
 492 {
 493     const   real sixth = 1.0/6.0;
 494     int          i, j;
 495     double       rsum, rav, sumaver, sumt;
 496
 497     sumaver = 0;
 498     sumt    = 0;
 499     for (i = 0; (i < nb); i++)
 500     {
 501         rsum = 0.0;
 502         rav  = 0.0;
 503         for (j = index[i]; (j < index[i+1]); j++)
 504         {
 505             if (viol)
 506             {
 507                 viol[j] += mypow(rt[j], -3.0);
 508             }
 509             rav     += sqr(rav3[j]);
 510             rsum    += mypow(rt[j], -6);
 511         }
 512         rsum    = std::max(0.0, mypow(rsum, -sixth)-bounds[i]);
 513         rav     = std::max(0.0, mypow(rav, -sixth)-bounds[i]);
 514
 515         sumt    += rsum;
 516         sumaver += rav;
 517     }
 518     *st = sumt;
 519     *sa = sumaver;
 520 }
 521
 522 static void analyse_disre(const char *voutfn,    int nframes,
 523                           real violaver[], real bounds[], int index[],
 524                           int pair[],      int nbounds,
 525                           const output_env_t oenv)
 526 {
 527     FILE   *vout;
 528     double  sum, sumt, sumaver;
 529     int     i, j;
 530
 531     /* Subtract bounds from distances, to calculate violations */
 532     calc_violations(violaver, violaver,
 533                     nbounds, pair, bounds, NULL, &sumt, &sumaver);
 534
 535 #ifdef DEBUG
 536     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 537             sumaver);
 538     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
 539             sumt);
 540 #endif
 541     vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
 542                     oenv);
 543     sum  = 0.0;
 544     sumt = 0.0;
 545     for (i = 0; (i < nbounds); i++)
 546     {
 547         /* Do ensemble averaging */
 548         sumaver = 0;
 549         for (j = pair[i]; (j < pair[i+1]); j++)
 550         {
 551             sumaver += sqr(violaver[j]/nframes);
 552         }
 553         sumaver = std::max(0.0, mypow(sumaver, minsixth)-bounds[i]);
 554
 555         sumt   += sumaver;
 556         sum     = std::max(sum, sumaver);
 557         fprintf(vout, "%10d  %10.5e\n", index[i], sumaver);
 558     }
 559 #ifdef DEBUG
 560     for (j = 0; (j < dr.ndr); j++)
 561     {
 562         fprintf(vout, "%10d  %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
 563     }
 564 #endif
 565     xvgrclose(vout);
 566
 567     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 568             sumt);
 569     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
 570
 571     do_view(oenv, voutfn, "-graphtype bar");
 572 }
 573
 574 static void einstein_visco(const char *fn, const char *fni, int nsets,
 575                            int nint, real **eneint,
 576                            real V, real T, double dt,
 577                            const output_env_t oenv)
 578 {
 579     FILE  *fp0, *fp1;
 580     double av[4], avold[4];
 581     double fac, di;
 582     int    i, j, m, nf4;
 583
 584     nf4 = nint/4 + 1;
 585
 586     for (i = 0; i <= nsets; i++)
 587     {
 588         avold[i] = 0;
 589     }
 590     fp0 = xvgropen(fni, "Shear viscosity integral",
 591                    "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
 592     fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
 593                    "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
 594     for (i = 0; i < nf4; i++)
 595     {
 596         for (m = 0; m <= nsets; m++)
 597         {
 598             av[m] = 0;
 599         }
 600         for (j = 0; j < nint - i; j++)
 601         {
 602             for (m = 0; m < nsets; m++)
 603             {
 604                 di   = sqr(eneint[m][j+i] - eneint[m][j]);
 605
 606                 av[m]     += di;
 607                 av[nsets] += di/nsets;
 608             }
 609         }
 610         /* Convert to SI for the viscosity */
 611         fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nint - i);
 612         fprintf(fp0, "%10g", i*dt);
 613         for (m = 0; (m <= nsets); m++)
 614         {
 615             av[m] = fac*av[m];
 616             fprintf(fp0, "  %10g", av[m]);
 617         }
 618         fprintf(fp0, "\n");
 619         fprintf(fp1, "%10g", (i + 0.5)*dt);
 620         for (m = 0; (m <= nsets); m++)
 621         {
 622             fprintf(fp1, "  %10g", (av[m]-avold[m])/dt);
 623             avold[m] = av[m];
 624         }
 625         fprintf(fp1, "\n");
 626     }
 627     xvgrclose(fp0);
 628     xvgrclose(fp1);
 629 }
 630
 631 typedef struct {
 632     gmx_int64_t     np;
 633     double          sum;
 634     double          sav;
 635     double          sav2;
 636 } ee_sum_t;
 637
 638 typedef struct {
 639     int             b;
 640     ee_sum_t        sum;
 641     gmx_int64_t     nst;
 642     gmx_int64_t     nst_min;
 643 } ener_ee_t;
 644
 645 static void clear_ee_sum(ee_sum_t *ees)
 646 {
 647     ees->sav  = 0;
 648     ees->sav2 = 0;
 649     ees->np   = 0;
 650     ees->sum  = 0;
 651 }
 652
 653 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
 654 {
 655     ees->np  += np;
 656     ees->sum += sum;
 657 }
 658
 659 static void add_ee_av(ee_sum_t *ees)
 660 {
 661     double av;
 662
 663     av         = ees->sum/ees->np;
 664     ees->sav  += av;
 665     ees->sav2 += av*av;
 666     ees->np    = 0;
 667     ees->sum   = 0;
 668 }
 669
 670 static double calc_ee2(int nb, ee_sum_t *ees)
 671 {
 672     return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
 673 }
 674
 675 static void set_ee_av(ener_ee_t *eee)
 676 {
 677     if (debug)
 678     {
 679         char buf[STEPSTRSIZE];
 680         fprintf(debug, "Storing average for err.est.: %s steps\n",
 681                 gmx_step_str(eee->nst, buf));
 682     }
 683     add_ee_av(&eee->sum);
 684     eee->b++;
 685     if (eee->b == 1 || eee->nst < eee->nst_min)
 686     {
 687         eee->nst_min = eee->nst;
 688     }
 689     eee->nst = 0;
 690 }
 691
 692 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
 693 {
 694     int             nb, i, f, nee;
 695     double          sum, sum2, sump, see2;
 696     gmx_int64_t     np, p, bound_nb;
 697     enerdat_t      *ed;
 698     exactsum_t     *es;
 699     gmx_bool        bAllZero;
 700     double          x, sx, sy, sxx, sxy;
 701     ener_ee_t      *eee;
 702
 703     /* Check if we have exact statistics over all points */
 704     for (i = 0; i < nset; i++)
 705     {
 706         ed             = &edat->s[i];
 707         ed->bExactStat = FALSE;
 708         if (edat->bHaveSums)
 709         {
 710             /* All energy file sum entries 0 signals no exact sums.
 711              * But if all energy values are 0, we still have exact sums.
 712              */
 713             bAllZero = TRUE;
 714             for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
 715             {
 716                 if (ed->ener[i] != 0)
 717                 {
 718                     bAllZero = FALSE;
 719                 }
 720                 ed->bExactStat = (ed->es[f].sum != 0);
 721             }
 722             if (bAllZero)
 723             {
 724                 ed->bExactStat = TRUE;
 725             }
 726         }
 727     }
 728
 729     snew(eee, nbmax+1);
 730     for (i = 0; i < nset; i++)
 731     {
 732         ed = &edat->s[i];
 733
 734         sum  = 0;
 735         sum2 = 0;
 736         np   = 0;
 737         sx   = 0;
 738         sy   = 0;
 739         sxx  = 0;
 740         sxy  = 0;
 741         for (nb = nbmin; nb <= nbmax; nb++)
 742         {
 743             eee[nb].b     = 0;
 744             clear_ee_sum(&eee[nb].sum);
 745             eee[nb].nst     = 0;
 746             eee[nb].nst_min = 0;
 747         }
 748         for (f = 0; f < edat->nframes; f++)
 749         {
 750             es = &ed->es[f];
 751
 752             if (ed->bExactStat)
 753             {
 754                 /* Add the sum and the sum of variances to the totals. */
 755                 p     = edat->points[f];
 756                 sump  = es->sum;
 757                 sum2 += es->sum2;
 758                 if (np > 0)
 759                 {
 760                     sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
 761                         *np*(np + p)/p;
 762                 }
 763             }
 764             else
 765             {
 766                 /* Add a single value to the sum and sum of squares. */
 767                 p     = 1;
 768                 sump  = ed->ener[f];
 769                 sum2 += dsqr(sump);
 770             }
 771
 772             /* sum has to be increased after sum2 */
 773             np  += p;
 774             sum += sump;
 775
 776             /* For the linear regression use variance 1/p.
 777              * Note that sump is the sum, not the average, so we don't need p*.
 778              */
 779             x    = edat->step[f] - 0.5*(edat->steps[f] - 1);
 780             sx  += p*x;
 781             sy  += sump;
 782             sxx += p*x*x;
 783             sxy += x*sump;
 784
 785             for (nb = nbmin; nb <= nbmax; nb++)
 786             {
 787                 /* Check if the current end step is closer to the desired
 788                  * block boundary than the next end step.
 789                  */
 790                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 791                 if (eee[nb].nst > 0 &&
 792                     bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
 793                 {
 794                     set_ee_av(&eee[nb]);
 795                 }
 796                 if (f == 0)
 797                 {
 798                     eee[nb].nst = 1;
 799                 }
 800                 else
 801                 {
 802                     eee[nb].nst += edat->step[f] - edat->step[f-1];
 803                 }
 804                 if (ed->bExactStat)
 805                 {
 806                     add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
 807                 }
 808                 else
 809                 {
 810                     add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
 811                 }
 812                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 813                 if (edat->step[f]*nb >= bound_nb)
 814                 {
 815                     set_ee_av(&eee[nb]);
 816                 }
 817             }
 818         }
 819
 820         edat->s[i].av = sum/np;
 821         if (ed->bExactStat)
 822         {
 823             edat->s[i].rmsd = std::sqrt(sum2/np);
 824         }
 825         else
 826         {
 827             edat->s[i].rmsd = std::sqrt(sum2/np - dsqr(edat->s[i].av));
 828         }
 829
 830         if (edat->nframes > 1)
 831         {
 832             edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
 833         }
 834         else
 835         {
 836             edat->s[i].slope = 0;
 837         }
 838
 839         nee  = 0;
 840         see2 = 0;
 841         for (nb = nbmin; nb <= nbmax; nb++)
 842         {
 843             /* Check if we actually got nb blocks and if the smallest
 844              * block is not shorter than 80% of the average.
 845              */
 846             if (debug)
 847             {
 848                 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
 849                 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
 850                         nb, eee[nb].b,
 851                         gmx_step_str(eee[nb].nst_min, buf1),
 852                         gmx_step_str(edat->nsteps, buf2));
 853             }
 854             if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
 855             {
 856                 see2 += calc_ee2(nb, &eee[nb].sum);
 857                 nee++;
 858             }
 859         }
 860         if (nee > 0)
 861         {
 862             edat->s[i].ee = std::sqrt(see2/nee);
 863         }
 864         else
 865         {
 866             edat->s[i].ee = -1;
 867         }
 868     }
 869     sfree(eee);
 870 }
 871
 872 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
 873 {
 874     enerdata_t *esum;
 875     enerdat_t  *s;
 876     int         f, i;
 877     double      sum;
 878
 879     snew(esum, 1);
 880     *esum = *edat;
 881     snew(esum->s, 1);
 882     s = &esum->s[0];
 883     snew(s->ener, esum->nframes);
 884     snew(s->es, esum->nframes);
 885
 886     s->bExactStat = TRUE;
 887     s->slope      = 0;
 888     for (i = 0; i < nset; i++)
 889     {
 890         if (!edat->s[i].bExactStat)
 891         {
 892             s->bExactStat = FALSE;
 893         }
 894         s->slope += edat->s[i].slope;
 895     }
 896
 897     for (f = 0; f < edat->nframes; f++)
 898     {
 899         sum = 0;
 900         for (i = 0; i < nset; i++)
 901         {
 902             sum += edat->s[i].ener[f];
 903         }
 904         s->ener[f] = sum;
 905         sum        = 0;
 906         for (i = 0; i < nset; i++)
 907         {
 908             sum += edat->s[i].es[f].sum;
 909         }
 910         s->es[f].sum  = sum;
 911         s->es[f].sum2 = 0;
 912     }
 913
 914     calc_averages(1, esum, nbmin, nbmax);
 915
 916     return esum;
 917 }
 918
 919 static char *ee_pr(double ee, char *buf)
 920 {
 921     char   tmp[100];
 922     double rnd;
 923
 924     if (ee < 0)
 925     {
 926         sprintf(buf, "%s", "--");
 927     }
 928     else
 929     {
 930         /* Round to two decimals by printing. */
 931         sprintf(tmp, "%.1e", ee);
 932         sscanf(tmp, "%lf", &rnd);
 933         sprintf(buf, "%g", rnd);
 934     }
 935
 936     return buf;
 937 }
 938
 939 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
 940 {
 941 /* Remove the drift by performing a fit to y = ax+b.
 942    Uses 5 iterations. */
 943     int    i, j, k;
 944     double delta;
 945
 946     edat->npoints = edat->nframes;
 947     edat->nsteps  = edat->nframes;
 948
 949     for (k = 0; (k < 5); k++)
 950     {
 951         for (i = 0; (i < nset); i++)
 952         {
 953             delta = edat->s[i].slope*dt;
 954
 955             if (NULL != debug)
 956             {
 957                 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
 958             }
 959
 960             for (j = 0; (j < edat->nframes); j++)
 961             {
 962                 edat->s[i].ener[j]   -= j*delta;
 963                 edat->s[i].es[j].sum  = 0;
 964                 edat->s[i].es[j].sum2 = 0;
 965             }
 966         }
 967         calc_averages(nset, edat, nbmin, nbmax);
 968     }
 969 }
 970
 971 static void calc_fluctuation_props(FILE *fp,
 972                                    gmx_bool bDriftCorr, real dt,
 973                                    int nset, int nmol,
 974                                    char **leg, enerdata_t *edat,
 975                                    int nbmin, int nbmax)
 976 {
 977     int    i, j;
 978     double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, varet;
 979     double NANO3;
 980     enum {
 981         eVol, eEnth, eTemp, eEtot, eNR
 982     };
 983     const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
 984     int         ii[eNR];
 985
 986     NANO3 = NANO*NANO*NANO;
 987     if (!bDriftCorr)
 988     {
 989         fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
 990                 "for spurious drift in the graphs. Note that this is not\n"
 991                 "a substitute for proper equilibration and sampling!\n");
 992     }
 993     else
 994     {
 995         remove_drift(nset, nbmin, nbmax, dt, edat);
 996     }
 997     for (i = 0; (i < eNR); i++)
 998     {
 999         for (ii[i] = 0; (ii[i] < nset &&
1000                          (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
1001         {
1002             ;
1003         }
1004 /*        if (ii[i] < nset)
1005             fprintf(fp,"Found %s data.\n",my_ener[i]);
1006  */ }
1007     /* Compute it all! */
1008     alpha = kappa = cp = dcp = cv = NOTSET;
1009
1010     /* Temperature */
1011     tt = NOTSET;
1012     if (ii[eTemp] < nset)
1013     {
1014         tt    = edat->s[ii[eTemp]].av;
1015     }
1016     /* Volume */
1017     vv = varv = NOTSET;
1018     if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1019     {
1020         vv    = edat->s[ii[eVol]].av*NANO3;
1021         varv  = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1022         kappa = (varv/vv)/(BOLTZMANN*tt);
1023     }
1024     /* Enthalpy */
1025     hh = varh = NOTSET;
1026     if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1027     {
1028         hh    = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1029         varh  = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1030         cp    = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1031     }
1032     /* Total energy */
1033     if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1034     {
1035         /* Only compute cv in constant volume runs, which we can test
1036            by checking whether the enthalpy was computed.
1037          */
1038         varet = sqr(edat->s[ii[eEtot]].rmsd);
1039         cv    = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1040     }
1041     /* Alpha, dcp */
1042     if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1043     {
1044         double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
1045         vh_sum = v_sum = h_sum = 0;
1046         for (j = 0; (j < edat->nframes); j++)
1047         {
1048             v       = edat->s[ii[eVol]].ener[j]*NANO3;
1049             h       = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
1050             v_sum  += v;
1051             h_sum  += h;
1052             vh_sum += (v*h);
1053         }
1054         vh_aver = vh_sum / edat->nframes;
1055         v_aver  = v_sum  / edat->nframes;
1056         h_aver  = h_sum  / edat->nframes;
1057         alpha   = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
1058         dcp     = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1059     }
1060
1061     if (tt != NOTSET)
1062     {
1063         if (nmol < 2)
1064         {
1065             fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1066                     nmol);
1067         }
1068         fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1069         fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1070         fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1071         fprintf(fp, "please use the g_dos program.\n\n");
1072         fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1073                 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1074
1075         if (debug != NULL)
1076         {
1077             if (varv != NOTSET)
1078             {
1079                 fprintf(fp, "varv  =  %10g (m^6)\n", varv*AVOGADRO/nmol);
1080             }
1081         }
1082         if (vv != NOTSET)
1083         {
1084             fprintf(fp, "Volume                                   = %10g m^3/mol\n",
1085                     vv*AVOGADRO/nmol);
1086         }
1087         if (varh != NOTSET)
1088         {
1089             fprintf(fp, "Enthalpy                                 = %10g kJ/mol\n",
1090                     hh*AVOGADRO/(KILO*nmol));
1091         }
1092         if (alpha != NOTSET)
1093         {
1094             fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
1095                     alpha);
1096         }
1097         if (kappa != NOTSET)
1098         {
1099             fprintf(fp, "Isothermal Compressibility Kappa         = %10g (J/m^3)\n",
1100                     kappa);
1101             fprintf(fp, "Adiabatic bulk modulus                   = %10g (m^3/J)\n",
1102                     1.0/kappa);
1103         }
1104         if (cp != NOTSET)
1105         {
1106             fprintf(fp, "Heat capacity at constant pressure Cp    = %10g J/mol K\n",
1107                     cp);
1108         }
1109         if (cv != NOTSET)
1110         {
1111             fprintf(fp, "Heat capacity at constant volume Cv      = %10g J/mol K\n",
1112                     cv);
1113         }
1114         if (dcp != NOTSET)
1115         {
1116             fprintf(fp, "Cp-Cv                                    =  %10g J/mol K\n",
1117                     dcp);
1118         }
1119         please_cite(fp, "Allen1987a");
1120     }
1121     else
1122     {
1123         fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
1124     }
1125 }
1126
1127 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
1128                          gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
1129                          gmx_bool bVisco, const char *visfn, int nmol,
1130                          gmx_int64_t start_step, double start_t,
1131                          gmx_int64_t step, double t,
1132                          real reftemp,
1133                          enerdata_t *edat,
1134                          int nset, int set[], gmx_bool *bIsEner,
1135                          char **leg, gmx_enxnm_t *enm,
1136                          real Vaver, real ezero,
1137                          int nbmin, int nbmax,
1138                          const output_env_t oenv)
1139 {
1140     FILE           *fp;
1141     /* Check out the printed manual for equations! */
1142     double          Dt, aver, stddev, errest, delta_t, totaldrift;
1143     enerdata_t     *esum = NULL;
1144     real            integral, intBulk, Temp = 0, Pres = 0;
1145     real            pr_aver, pr_stddev, pr_errest;
1146     double          beta = 0, expE, expEtot, *fee = NULL;
1147     gmx_int64_t     nsteps;
1148     int             nexact, nnotexact;
1149     int             i, j, nout;
1150     char            buf[256], eebuf[100];
1151
1152     nsteps  = step - start_step + 1;
1153     if (nsteps < 1)
1154     {
1155         fprintf(stdout, "Not enough steps (%s) for statistics\n",
1156                 gmx_step_str(nsteps, buf));
1157     }
1158     else
1159     {
1160         /* Calculate the time difference */
1161         delta_t = t - start_t;
1162
1163         fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1164                 gmx_step_str(nsteps, buf), start_t, t, nset);
1165
1166         calc_averages(nset, edat, nbmin, nbmax);
1167
1168         if (bSum)
1169         {
1170             esum = calc_sum(nset, edat, nbmin, nbmax);
1171         }
1172
1173         if (!edat->bHaveSums)
1174         {
1175             nexact    = 0;
1176             nnotexact = nset;
1177         }
1178         else
1179         {
1180             nexact    = 0;
1181             nnotexact = 0;
1182             for (i = 0; (i < nset); i++)
1183             {
1184                 if (edat->s[i].bExactStat)
1185                 {
1186                     nexact++;
1187                 }
1188                 else
1189                 {
1190                     nnotexact++;
1191                 }
1192             }
1193         }
1194
1195         if (nnotexact == 0)
1196         {
1197             fprintf(stdout, "All statistics are over %s points\n",
1198                     gmx_step_str(edat->npoints, buf));
1199         }
1200         else if (nexact == 0 || edat->npoints == edat->nframes)
1201         {
1202             fprintf(stdout, "All statistics are over %d points (frames)\n",
1203                     edat->nframes);
1204         }
1205         else
1206         {
1207             fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
1208             for (i = 0; (i < nset); i++)
1209             {
1210                 if (!edat->s[i].bExactStat)
1211                 {
1212                     fprintf(stdout, " '%s'", leg[i]);
1213                 }
1214             }
1215             fprintf(stdout, " %s has statistics over %d points (frames)\n",
1216                     nnotexact == 1 ? "is" : "are", edat->nframes);
1217             fprintf(stdout, "All other statistics are over %s points\n",
1218                     gmx_step_str(edat->npoints, buf));
1219         }
1220         fprintf(stdout, "\n");
1221
1222         fprintf(stdout, "%-24s %10s %10s %10s %10s",
1223                 "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
1224         if (bFee)
1225         {
1226             fprintf(stdout, "  %10s\n", "-kT ln<e^(E/kT)>");
1227         }
1228         else
1229         {
1230             fprintf(stdout, "\n");
1231         }
1232         fprintf(stdout, "-------------------------------------------------------------------------------\n");
1233
1234         /* Initiate locals, only used with -sum */
1235         expEtot = 0;
1236         if (bFee)
1237         {
1238             beta = 1.0/(BOLTZ*reftemp);
1239             snew(fee, nset);
1240         }
1241         for (i = 0; (i < nset); i++)
1242         {
1243             aver   = edat->s[i].av;
1244             stddev = edat->s[i].rmsd;
1245             errest = edat->s[i].ee;
1246
1247             if (bFee)
1248             {
1249                 expE = 0;
1250                 for (j = 0; (j < edat->nframes); j++)
1251                 {
1252                     expE += std::exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1253                 }
1254                 if (bSum)
1255                 {
1256                     expEtot += expE/edat->nframes;
1257                 }
1258
1259                 fee[i] = std::log(expE/edat->nframes)/beta + aver/nmol;
1260             }
1261             if (std::strstr(leg[i], "empera") != NULL)
1262             {
1263                 Temp = aver;
1264             }
1265             else if (std::strstr(leg[i], "olum") != NULL)
1266             {
1267                 Vaver = aver;
1268             }
1269             else if (std::strstr(leg[i], "essure") != NULL)
1270             {
1271                 Pres = aver;
1272             }
1273             if (bIsEner[i])
1274             {
1275                 pr_aver   = aver/nmol-ezero;
1276                 pr_stddev = stddev/nmol;
1277                 pr_errest = errest/nmol;
1278             }
1279             else
1280             {
1281                 pr_aver   = aver;
1282                 pr_stddev = stddev;
1283                 pr_errest = errest;
1284             }
1285
1286             /* Multiply the slope in steps with the number of steps taken */
1287             totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1288             if (bIsEner[i])
1289             {
1290                 totaldrift /= nmol;
1291             }
1292
1293             fprintf(stdout, "%-24s %10g %10s %10g %10g",
1294                     leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift);
1295             if (bFee)
1296             {
1297                 fprintf(stdout, "  %10g", fee[i]);
1298             }
1299
1300             fprintf(stdout, "  (%s)\n", enm[set[i]].unit);
1301
1302             if (bFluct)
1303             {
1304                 for (j = 0; (j < edat->nframes); j++)
1305                 {
1306                     edat->s[i].ener[j] -= aver;
1307                 }
1308             }
1309         }
1310         if (bSum)
1311         {
1312             totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1313             fprintf(stdout, "%-24s %10g %10s %10s %10g  (%s)",
1314                     "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf),
1315                     "--", totaldrift/nmol, enm[set[0]].unit);
1316             /* pr_aver,pr_stddev,a,totaldrift */
1317             if (bFee)
1318             {
1319                 fprintf(stdout, "  %10g  %10g\n",
1320                         std::log(expEtot)/beta + esum->s[0].av/nmol, std::log(expEtot)/beta);
1321             }
1322             else
1323             {
1324                 fprintf(stdout, "\n");
1325             }
1326         }
1327
1328         /* Do correlation function */
1329         if (edat->nframes > 1)
1330         {
1331             Dt = delta_t/(edat->nframes - 1);
1332         }
1333         else
1334         {
1335             Dt = 0;
1336         }
1337         if (bVisco)
1338         {
1339             const char* leg[] = { "Shear", "Bulk" };
1340             real        factor;
1341             real      **eneset;
1342             real      **eneint;
1343
1344             /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1345
1346             /* Symmetrise tensor! (and store in first three elements)
1347              * And subtract average pressure!
1348              */
1349             snew(eneset, 12);
1350             for (i = 0; i < 12; i++)
1351             {
1352                 snew(eneset[i], edat->nframes);
1353             }
1354             for (i = 0; (i < edat->nframes); i++)
1355             {
1356                 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1357                 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1358                 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1359                 for (j = 3; j <= 11; j++)
1360                 {
1361                     eneset[j][i] = edat->s[j].ener[i];
1362                 }
1363                 eneset[11][i] -= Pres;
1364             }
1365
1366             /* Determine integrals of the off-diagonal pressure elements */
1367             snew(eneint, 3);
1368             for (i = 0; i < 3; i++)
1369             {
1370                 snew(eneint[i], edat->nframes + 1);
1371             }
1372             eneint[0][0] = 0;
1373             eneint[1][0] = 0;
1374             eneint[2][0] = 0;
1375             for (i = 0; i < edat->nframes; i++)
1376             {
1377                 eneint[0][i+1] = eneint[0][i] + 0.5*(edat->s[1].es[i].sum + edat->s[3].es[i].sum)*Dt/edat->points[i];
1378                 eneint[1][i+1] = eneint[1][i] + 0.5*(edat->s[2].es[i].sum + edat->s[6].es[i].sum)*Dt/edat->points[i];
1379                 eneint[2][i+1] = eneint[2][i] + 0.5*(edat->s[5].es[i].sum + edat->s[7].es[i].sum)*Dt/edat->points[i];
1380             }
1381
1382             einstein_visco("evisco.xvg", "eviscoi.xvg",
1383                            3, edat->nframes+1, eneint, Vaver, Temp, Dt, oenv);
1384
1385             for (i = 0; i < 3; i++)
1386             {
1387                 sfree(eneint[i]);
1388             }
1389             sfree(eneint);
1390
1391             /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1392             /* Do it for shear viscosity */
1393             std::strcpy(buf, "Shear Viscosity");
1394             low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
1395                             (edat->nframes+1)/2, eneset, Dt,
1396                             eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1397
1398             /* Now for bulk viscosity */
1399             std::strcpy(buf, "Bulk Viscosity");
1400             low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1401                             (edat->nframes+1)/2, &(eneset[11]), Dt,
1402                             eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1403
1404             factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1405             fp     = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1406             xvgr_legend(fp, asize(leg), leg, oenv);
1407
1408             /* Use trapezium rule for integration */
1409             integral = 0;
1410             intBulk  = 0;
1411             nout     = get_acfnout();
1412             if ((nout < 2) || (nout >= edat->nframes/2))
1413             {
1414                 nout = edat->nframes/2;
1415             }
1416             for (i = 1; (i < nout); i++)
1417             {
1418                 integral += 0.5*(eneset[0][i-1]  + eneset[0][i])*factor;
1419                 intBulk  += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1420                 fprintf(fp, "%10g  %10g  %10g\n", (i*Dt), integral, intBulk);
1421             }
1422             xvgrclose(fp);
1423
1424             for (i = 0; i < 12; i++)
1425             {
1426                 sfree(eneset[i]);
1427             }
1428             sfree(eneset);
1429         }
1430         else if (bCorr)
1431         {
1432             if (bFluct)
1433             {
1434                 std::strcpy(buf, "Autocorrelation of Energy Fluctuations");
1435             }
1436             else
1437             {
1438                 std::strcpy(buf, "Energy Autocorrelation");
1439             }
1440 #if 0
1441             do_autocorr(corrfn, oenv, buf, edat->nframes,
1442                         bSum ? 1                 : nset,
1443                         bSum ? &edat->s[nset-1].ener : eneset,
1444                         (delta_t/edat->nframes), eacNormal, FALSE);
1445 #endif
1446         }
1447     }
1448 }
1449
1450 static void print_time(FILE *fp, double t)
1451 {
1452     fprintf(fp, "%12.6f", t);
1453 }
1454
1455 static void print1(FILE *fp, gmx_bool bDp, real e)
1456 {
1457     if (bDp)
1458     {
1459         fprintf(fp, "  %16.12f", e);
1460     }
1461     else
1462     {
1463         fprintf(fp, "  %10.6f", e);
1464     }
1465 }
1466
1467 static void fec(const char *ene2fn, const char *runavgfn,
1468                 real reftemp, int nset, int set[], char *leg[],
1469                 enerdata_t *edat, double time[],
1470                 const output_env_t oenv)
1471 {
1472     const char * ravgleg[] = {
1473         "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1474         "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
1475     };
1476     FILE        *fp;
1477     ener_file_t  enx;
1478     int          timecheck, nenergy, nenergy2, maxenergy;
1479     int          i, j;
1480     gmx_bool     bCont;
1481     real         aver, beta;
1482     real       **eneset2;
1483     double       dE, sum;
1484     gmx_enxnm_t *enm = NULL;
1485     t_enxframe  *fr;
1486     char         buf[22];
1487
1488     /* read second energy file */
1489     snew(fr, 1);
1490     enm = NULL;
1491     enx = open_enx(ene2fn, "r");
1492     do_enxnms(enx, &(fr->nre), &enm);
1493
1494     snew(eneset2, nset+1);
1495     nenergy2  = 0;
1496     maxenergy = 0;
1497     timecheck = 0;
1498     do
1499     {
1500         /* This loop searches for the first frame (when -b option is given),
1501          * or when this has been found it reads just one energy frame
1502          */
1503         do
1504         {
1505             bCont = do_enx(enx, fr);
1506
1507             if (bCont)
1508             {
1509                 timecheck = check_times(fr->t);
1510             }
1511
1512         }
1513         while (bCont && (timecheck < 0));
1514
1515         /* Store energies for analysis afterwards... */
1516         if ((timecheck == 0) && bCont)
1517         {
1518             if (fr->nre > 0)
1519             {
1520                 if (nenergy2 >= maxenergy)
1521                 {
1522                     maxenergy += 1000;
1523                     for (i = 0; i <= nset; i++)
1524                     {
1525                         srenew(eneset2[i], maxenergy);
1526                     }
1527                 }
1528                 GMX_RELEASE_ASSERT(time != NULL, "trying to dereference NULL time pointer");
1529
1530                 if (fr->t != time[nenergy2])
1531                 {
1532                     fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
1533                             fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
1534                 }
1535                 for (i = 0; i < nset; i++)
1536                 {
1537                     eneset2[i][nenergy2] = fr->ener[set[i]].e;
1538                 }
1539                 nenergy2++;
1540             }
1541         }
1542     }
1543     while (bCont && (timecheck == 0));
1544
1545     /* check */
1546     if (edat->nframes != nenergy2)
1547     {
1548         fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
1549                 edat->nframes, nenergy2);
1550     }
1551     nenergy = std::min(edat->nframes, nenergy2);
1552
1553     /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1554     fp = NULL;
1555     if (runavgfn)
1556     {
1557         fp = xvgropen(runavgfn, "Running average free energy difference",
1558                       "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
1559         xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1560     }
1561     fprintf(stdout, "\n%-24s %10s\n",
1562             "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1563     sum  = 0;
1564     beta = 1.0/(BOLTZ*reftemp);
1565     for (i = 0; i < nset; i++)
1566     {
1567         if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1568         {
1569             fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
1570                     leg[i], enm[set[i]].name);
1571         }
1572         for (j = 0; j < nenergy; j++)
1573         {
1574             dE   = eneset2[i][j] - edat->s[i].ener[j];
1575             sum += std::exp(-dE*beta);
1576             if (fp)
1577             {
1578                 fprintf(fp, "%10g %10g %10g\n",
1579                         time[j], dE, -BOLTZ*reftemp*std::log(sum/(j+1)) );
1580             }
1581         }
1582         aver = -BOLTZ*reftemp*std::log(sum/nenergy);
1583         fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1584     }
1585     if (fp)
1586     {
1587         xvgrclose(fp);
1588     }
1589     sfree(fr);
1590 }
1591
1592
1593 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
1594                     const char *filename, gmx_bool bDp,
1595                     int *blocks, int *hists, int *samples, int *nlambdas,
1596                     const output_env_t oenv)
1597 {
1598     const char  *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1599     char         title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1600     char         buf[STRLEN];
1601     int          nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1602     int          i, j, k;
1603     /* coll data */
1604     double       temp              = 0, start_time = 0, delta_time = 0, start_lambda = 0;
1605     static int   setnr             = 0;
1606     double      *native_lambda_vec = NULL;
1607     const char **lambda_components = NULL;
1608     int          n_lambda_vec      = 0;
1609
1610     /* now count the blocks & handle the global dh data */
1611     for (i = 0; i < fr->nblock; i++)
1612     {
1613         if (fr->block[i].id == enxDHHIST)
1614         {
1615             nblock_hist++;
1616         }
1617         else if (fr->block[i].id == enxDH)
1618         {
1619             nblock_dh++;
1620         }
1621         else if (fr->block[i].id == enxDHCOLL)
1622         {
1623             nblock_dhcoll++;
1624             if ( (fr->block[i].nsub < 1) ||
1625                  (fr->block[i].sub[0].type != xdr_datatype_double) ||
1626                  (fr->block[i].sub[0].nr < 5))
1627             {
1628                 gmx_fatal(FARGS, "Unexpected block data");
1629             }
1630
1631             /* read the data from the DHCOLL block */
1632             temp            = fr->block[i].sub[0].dval[0];
1633             start_time      = fr->block[i].sub[0].dval[1];
1634             delta_time      = fr->block[i].sub[0].dval[2];
1635             start_lambda    = fr->block[i].sub[0].dval[3];
1636             if (fr->block[i].nsub > 1)
1637             {
1638                 if (n_lambda_vec == 0)
1639                 {
1640                     n_lambda_vec = fr->block[i].sub[1].ival[1];
1641                 }
1642                 else
1643                 {
1644                     if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1645                     {
1646                         gmx_fatal(FARGS,
1647                                   "Unexpected change of basis set in lambda");
1648                     }
1649                 }
1650                 if (lambda_components == NULL)
1651                 {
1652                     snew(lambda_components, n_lambda_vec);
1653                 }
1654                 if (native_lambda_vec == NULL)
1655                 {
1656                     snew(native_lambda_vec, n_lambda_vec);
1657                 }
1658                 for (j = 0; j < n_lambda_vec; j++)
1659                 {
1660                     native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1661                     lambda_components[j] =
1662                         efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1663                 }
1664             }
1665         }
1666     }
1667
1668     if (nblock_hist == 0 && nblock_dh == 0)
1669     {
1670         /* don't do anything */
1671         return;
1672     }
1673     if (nblock_hist > 0 && nblock_dh > 0)
1674     {
1675         gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1676     }
1677     if (!*fp_dhdl)
1678     {
1679         if (nblock_dh > 0)
1680         {
1681             /* we have standard, non-histogram data --
1682                call open_dhdl to open the file */
1683             /* TODO this is an ugly hack that needs to be fixed: this will only
1684                work if the order of data is always the same and if we're
1685                only using the g_energy compiled with the mdrun that produced
1686                the ener.edr. */
1687             *fp_dhdl = open_dhdl(filename, ir, oenv);
1688         }
1689         else
1690         {
1691             sprintf(title, "N(%s)", deltag);
1692             sprintf(label_x, "%s (%s)", deltag, unit_energy);
1693             sprintf(label_y, "Samples");
1694             *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1695             sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1696             xvgr_subtitle(*fp_dhdl, buf, oenv);
1697         }
1698     }
1699
1700     (*hists)   += nblock_hist;
1701     (*blocks)  += nblock_dh;
1702     (*nlambdas) = nblock_hist+nblock_dh;
1703
1704     /* write the data */
1705     if (nblock_hist > 0)
1706     {
1707         gmx_int64_t sum = 0;
1708         /* histograms */
1709         for (i = 0; i < fr->nblock; i++)
1710         {
1711             t_enxblock *blk = &(fr->block[i]);
1712             if (blk->id == enxDHHIST)
1713             {
1714                 double          foreign_lambda, dx;
1715                 gmx_int64_t     x0;
1716                 int             nhist, derivative;
1717
1718                 /* check the block types etc. */
1719                 if ( (blk->nsub < 2) ||
1720                      (blk->sub[0].type != xdr_datatype_double) ||
1721                      (blk->sub[1].type != xdr_datatype_int64) ||
1722                      (blk->sub[0].nr < 2)  ||
1723                      (blk->sub[1].nr < 2) )
1724                 {
1725                     gmx_fatal(FARGS, "Unexpected block data in file");
1726                 }
1727                 foreign_lambda = blk->sub[0].dval[0];
1728                 dx             = blk->sub[0].dval[1];
1729                 nhist          = blk->sub[1].lval[0];
1730                 derivative     = blk->sub[1].lval[1];
1731                 for (j = 0; j < nhist; j++)
1732                 {
1733                     const char *lg[1];
1734                     x0 = blk->sub[1].lval[2+j];
1735
1736                     if (!derivative)
1737                     {
1738                         sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1739                                 deltag, lambda, foreign_lambda,
1740                                 lambda, start_lambda);
1741                     }
1742                     else
1743                     {
1744                         sprintf(legend, "N(%s | %s=%g)",
1745                                 dhdl, lambda, start_lambda);
1746                     }
1747
1748                     lg[0] = legend;
1749                     xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1750                     setnr++;
1751                     for (k = 0; k < blk->sub[j+2].nr; k++)
1752                     {
1753                         int    hist;
1754                         double xmin, xmax;
1755
1756                         hist = blk->sub[j+2].ival[k];
1757                         xmin = (x0+k)*dx;
1758                         xmax = (x0+k+1)*dx;
1759                         fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1760                                 xmax, hist);
1761                         sum += hist;
1762                     }
1763                     /* multiple histogram data blocks in one histogram
1764                        mean that the second one is the reverse of the first one:
1765                        for dhdl derivatives, it's important to know both the
1766                        maximum and minimum values */
1767                     dx = -dx;
1768                 }
1769             }
1770         }
1771         (*samples) += static_cast<int>(sum/nblock_hist);
1772     }
1773     else
1774     {
1775         /* raw dh */
1776         int    len      = 0;
1777
1778         for (i = 0; i < fr->nblock; i++)
1779         {
1780             t_enxblock *blk = &(fr->block[i]);
1781             if (blk->id == enxDH)
1782             {
1783                 if (len == 0)
1784                 {
1785                     len = blk->sub[2].nr;
1786                 }
1787                 else
1788                 {
1789                     if (len != blk->sub[2].nr)
1790                     {
1791                         gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1792                     }
1793                 }
1794             }
1795         }
1796         (*samples) += len;
1797
1798         for (i = 0; i < len; i++)
1799         {
1800             double time = start_time + delta_time*i;
1801
1802             fprintf(*fp_dhdl, "%.4f ", time);
1803
1804             for (j = 0; j < fr->nblock; j++)
1805             {
1806                 t_enxblock *blk = &(fr->block[j]);
1807                 if (blk->id == enxDH)
1808                 {
1809                     double value;
1810                     if (blk->sub[2].type == xdr_datatype_float)
1811                     {
1812                         value = blk->sub[2].fval[i];
1813                     }
1814                     else
1815                     {
1816                         value = blk->sub[2].dval[i];
1817                     }
1818                     /* we need to decide which data type it is based on the count*/
1819
1820                     if (j == 1 && ir->bExpanded)
1821                     {
1822                         fprintf(*fp_dhdl, "%4d", static_cast<int>(value));   /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
1823                     }
1824                     else
1825                     {
1826                         if (bDp)
1827                         {
1828                             fprintf(*fp_dhdl, " %#.12g", value);   /* print normal precision */
1829                         }
1830                         else
1831                         {
1832                             fprintf(*fp_dhdl, " %#.8g", value);   /* print normal precision */
1833                         }
1834                     }
1835                 }
1836             }
1837             fprintf(*fp_dhdl, "\n");
1838         }
1839     }
1840 }
1841
1842
1843 int gmx_energy(int argc, char *argv[])
1844 {
1845     const char        *desc[] = {
1846         "[THISMODULE] extracts energy components or distance restraint",
1847         "data from an energy file. The user is prompted to interactively",
1848         "select the desired energy terms.[PAR]",
1849
1850         "Average, RMSD, and drift are calculated with full precision from the",
1851         "simulation (see printed manual). Drift is calculated by performing",
1852         "a least-squares fit of the data to a straight line. The reported total drift",
1853         "is the difference of the fit at the first and last point.",
1854         "An error estimate of the average is given based on a block averages",
1855         "over 5 blocks using the full-precision averages. The error estimate",
1856         "can be performed over multiple block lengths with the options",
1857         "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1858         "[BB]Note[bb] that in most cases the energy files contains averages over all",
1859         "MD steps, or over many more points than the number of frames in",
1860         "energy file. This makes the [THISMODULE] statistics output more accurate",
1861         "than the [REF].xvg[ref] output. When exact averages are not present in the energy",
1862         "file, the statistics mentioned above are simply over the single, per-frame",
1863         "energy values.[PAR]",
1864
1865         "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1866
1867         "Some fluctuation-dependent properties can be calculated provided",
1868         "the correct energy terms are selected, and that the command line option",
1869         "[TT]-fluct_props[tt] is given. The following properties",
1870         "will be computed:",
1871         "",
1872         "===============================  ===================",
1873         "Property                         Energy terms needed",
1874         "===============================  ===================",
1875         "Heat capacity C[SUB]p[sub] (NPT sims):    Enthalpy, Temp",
1876         "Heat capacity C[SUB]v[sub] (NVT sims):    Etot, Temp",
1877         "Thermal expansion coeff. (NPT):  Enthalpy, Vol, Temp",
1878         "Isothermal compressibility:      Vol, Temp",
1879         "Adiabatic bulk modulus:          Vol, Temp",
1880         "===============================  ===================",
1881         "",
1882         "You always need to set the number of molecules [TT]-nmol[tt].",
1883         "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1884         "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]"
1885         "When the [TT]-viol[tt] option is set, the time averaged",
1886         "violations are plotted and the running time-averaged and",
1887         "instantaneous sum of violations are recalculated. Additionally",
1888         "running time-averaged and instantaneous distances between",
1889         "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1890
1891         "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1892         "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1893         "The first two options plot the orientation, the last three the",
1894         "deviations of the orientations from the experimental values.",
1895         "The options that end on an 'a' plot the average over time",
1896         "as a function of restraint. The options that end on a 't'",
1897         "prompt the user for restraint label numbers and plot the data",
1898         "as a function of time. Option [TT]-odr[tt] plots the RMS",
1899         "deviation as a function of restraint.",
1900         "When the run used time or ensemble averaged orientation restraints,",
1901         "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1902         "not ensemble-averaged orientations and deviations instead of",
1903         "the time and ensemble averages.[PAR]",
1904
1905         "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1906         "tensor for each orientation restraint experiment. With option",
1907         "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1908
1909         "Option [TT]-odh[tt] extracts and plots the free energy data",
1910         "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1911         "from the [TT]ener.edr[tt] file.[PAR]",
1912
1913         "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1914         "difference with an ideal gas state::",
1915         "",
1916         "  [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1917         "  [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1918         "",
1919         "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1920         "the average is over the ensemble (or time in a trajectory).",
1921         "Note that this is in principle",
1922         "only correct when averaging over the whole (Boltzmann) ensemble",
1923         "and using the potential energy. This also allows for an entropy",
1924         "estimate using::",
1925         "",
1926         "  [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T",
1927         "  [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
1928         "",
1929
1930         "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1931         "difference is calculated::",
1932         "",
1933         "  dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1934         "",
1935         "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1936         "files, and the average is over the ensemble A. The running average",
1937         "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1938         "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1939
1940     };
1941     static gmx_bool    bSum    = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1942     static gmx_bool    bDp     = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1943     static int         skip    = 0, nmol = 1, nbmin = 5, nbmax = 5;
1944     static real        reftemp = 300.0, ezero = 0;
1945     t_pargs            pa[]    = {
1946         { "-fee",   FALSE, etBOOL,  {&bFee},
1947           "Do a free energy estimate" },
1948         { "-fetemp", FALSE, etREAL, {&reftemp},
1949           "Reference temperature for free energy calculation" },
1950         { "-zero", FALSE, etREAL, {&ezero},
1951           "Subtract a zero-point energy" },
1952         { "-sum",  FALSE, etBOOL, {&bSum},
1953           "Sum the energy terms selected rather than display them all" },
1954         { "-dp",   FALSE, etBOOL, {&bDp},
1955           "Print energies in high precision" },
1956         { "-nbmin", FALSE, etINT, {&nbmin},
1957           "Minimum number of blocks for error estimate" },
1958         { "-nbmax", FALSE, etINT, {&nbmax},
1959           "Maximum number of blocks for error estimate" },
1960         { "-mutot", FALSE, etBOOL, {&bMutot},
1961           "Compute the total dipole moment from the components" },
1962         { "-skip", FALSE, etINT,  {&skip},
1963           "Skip number of frames between data points" },
1964         { "-aver", FALSE, etBOOL, {&bPrAll},
1965           "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1966         { "-nmol", FALSE, etINT,  {&nmol},
1967           "Number of molecules in your sample: the energies are divided by this number" },
1968         { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1969           "Compute properties based on energy fluctuations, like heat capacity" },
1970         { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1971           "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1972         { "-fluc", FALSE, etBOOL, {&bFluct},
1973           "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1974         { "-orinst", FALSE, etBOOL, {&bOrinst},
1975           "Analyse instantaneous orientation data" },
1976         { "-ovec", FALSE, etBOOL, {&bOvec},
1977           "Also plot the eigenvectors with [TT]-oten[tt]" }
1978     };
1979     const char       * drleg[] = {
1980         "Running average",
1981         "Instantaneous"
1982     };
1983     static const char *setnm[] = {
1984         "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1985         "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1986         "Volume",  "Pressure"
1987     };
1988
1989     FILE              *out     = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1990     FILE              *fp_dhdl = NULL;
1991     ener_file_t        fp;
1992     int                timecheck = 0;
1993     gmx_mtop_t         mtop;
1994     gmx_localtop_t    *top = NULL;
1995     t_inputrec         ir;
1996     enerdata_t         edat;
1997     gmx_enxnm_t       *enm = NULL;
1998     t_enxframe        *frame, *fr = NULL;
1999     int                cur = 0;
2000 #define NEXT (1-cur)
2001     int                nre, teller, teller_disre, nfr;
2002     gmx_int64_t        start_step;
2003     int                nor = 0, nex = 0, norfr = 0, enx_i = 0;
2004     real               start_t;
2005     real              *bounds  = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
2006     int               *index   = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
2007     int                nbounds = 0, npairs;
2008     gmx_bool           bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
2009     gmx_bool           bFoundStart, bCont, bVisco;
2010     double             sum, sumaver, sumt, dbl;
2011     double            *time = NULL;
2012     real               Vaver;
2013     int               *set     = NULL, i, j, k, nset, sss;
2014     gmx_bool          *bIsEner = NULL;
2015     char             **pairleg, **odtleg, **otenleg;
2016     char             **leg = NULL;
2017     char              *anm_j, *anm_k, *resnm_j, *resnm_k;
2018     int                resnr_j, resnr_k;
2019     const char        *orinst_sub = "@ subtitle \"instantaneous\"\n";
2020     char               buf[256];
2021     output_env_t       oenv;
2022     t_enxblock        *blk       = NULL;
2023     t_enxblock        *blk_disre = NULL;
2024     int                ndisre    = 0;
2025     int                dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2026
2027     t_filenm           fnm[] = {
2028         { efEDR, "-f",    NULL,      ffREAD  },
2029         { efEDR, "-f2",   NULL,      ffOPTRD },
2030         { efTPR, "-s",    NULL,      ffOPTRD },
2031         { efXVG, "-o",    "energy",  ffWRITE },
2032         { efXVG, "-viol", "violaver", ffOPTWR },
2033         { efXVG, "-pairs", "pairs",   ffOPTWR },
2034         { efXVG, "-ora",  "orienta", ffOPTWR },
2035         { efXVG, "-ort",  "orientt", ffOPTWR },
2036         { efXVG, "-oda",  "orideva", ffOPTWR },
2037         { efXVG, "-odr",  "oridevr", ffOPTWR },
2038         { efXVG, "-odt",  "oridevt", ffOPTWR },
2039         { efXVG, "-oten", "oriten",  ffOPTWR },
2040         { efXVG, "-corr", "enecorr", ffOPTWR },
2041         { efXVG, "-vis",  "visco",   ffOPTWR },
2042         { efXVG, "-ravg", "runavgdf", ffOPTWR },
2043         { efXVG, "-odh",  "dhdl", ffOPTWR }
2044     };
2045 #define NFILE asize(fnm)
2046     int                npargs;
2047     t_pargs           *ppa;
2048
2049     npargs = asize(pa);
2050     ppa    = add_acf_pargs(&npargs, pa);
2051     if (!parse_common_args(&argc, argv,
2052                            PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END,
2053                            NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
2054     {
2055         return 0;
2056     }
2057
2058     bDRAll = opt2bSet("-pairs", NFILE, fnm);
2059     bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2060     bORA   = opt2bSet("-ora", NFILE, fnm);
2061     bORT   = opt2bSet("-ort", NFILE, fnm);
2062     bODA   = opt2bSet("-oda", NFILE, fnm);
2063     bODR   = opt2bSet("-odr", NFILE, fnm);
2064     bODT   = opt2bSet("-odt", NFILE, fnm);
2065     bORIRE = bORA || bORT || bODA || bODR || bODT;
2066     bOTEN  = opt2bSet("-oten", NFILE, fnm);
2067     bDHDL  = opt2bSet("-odh", NFILE, fnm);
2068
2069     nset = 0;
2070
2071     snew(frame, 2);
2072     fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2073     do_enxnms(fp, &nre, &enm);
2074
2075     Vaver = -1;
2076
2077     bVisco = opt2bSet("-vis", NFILE, fnm);
2078
2079     if ((!bDisRe) && (!bDHDL))
2080     {
2081         if (bVisco)
2082         {
2083             nset = asize(setnm);
2084             snew(set, nset);
2085             /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2086             for (j = 0; j < nset; j++)
2087             {
2088                 for (i = 0; i < nre; i++)
2089                 {
2090                     if (std::strstr(enm[i].name, setnm[j]))
2091                     {
2092                         set[j] = i;
2093                         break;
2094                     }
2095                 }
2096                 if (i == nre)
2097                 {
2098                     if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2099                     {
2100                         printf("Enter the box volume (" unit_volume "): ");
2101                         if (1 != scanf("%lf", &dbl))
2102                         {
2103                             gmx_fatal(FARGS, "Error reading user input");
2104                         }
2105                         Vaver = dbl;
2106                     }
2107                     else
2108                     {
2109                         gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2110                                   setnm[j]);
2111                     }
2112                 }
2113             }
2114         }
2115         else
2116         {
2117             set = select_by_name(nre, enm, &nset);
2118         }
2119         /* Print all the different units once */
2120         sprintf(buf, "(%s)", enm[set[0]].unit);
2121         for (i = 1; i < nset; i++)
2122         {
2123             for (j = 0; j < i; j++)
2124             {
2125                 if (std::strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2126                 {
2127                     break;
2128                 }
2129             }
2130             if (j == i)
2131             {
2132                 std::strcat(buf, ", (");
2133                 std::strcat(buf, enm[set[i]].unit);
2134                 std::strcat(buf, ")");
2135             }
2136         }
2137         out = xvgropen(opt2fn("-o", NFILE, fnm), "GROMACS Energies", "Time (ps)", buf,
2138                        oenv);
2139
2140         snew(leg, nset+1);
2141         for (i = 0; (i < nset); i++)
2142         {
2143             leg[i] = enm[set[i]].name;
2144         }
2145         if (bSum)
2146         {
2147             leg[nset] = gmx_strdup("Sum");
2148             xvgr_legend(out, nset+1, (const char**)leg, oenv);
2149         }
2150         else
2151         {
2152             xvgr_legend(out, nset, (const char**)leg, oenv);
2153         }
2154
2155         snew(bIsEner, nset);
2156         for (i = 0; (i < nset); i++)
2157         {
2158             bIsEner[i] = FALSE;
2159             for (j = 0; (j <= F_ETOT); j++)
2160             {
2161                 bIsEner[i] = bIsEner[i] ||
2162                     (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2163             }
2164         }
2165
2166         if (bPrAll && nset > 1)
2167         {
2168             gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2169         }
2170
2171         time = NULL;
2172
2173         if (bORIRE || bOTEN)
2174         {
2175             get_orires_parms(ftp2fn(efTPR, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2176         }
2177
2178         if (bORIRE)
2179         {
2180             if (bOrinst)
2181             {
2182                 enx_i = enxORI;
2183             }
2184             else
2185             {
2186                 enx_i = enxOR;
2187             }
2188
2189             if (bORA || bODA)
2190             {
2191                 snew(orient, nor);
2192             }
2193             if (bODR)
2194             {
2195                 snew(odrms, nor);
2196             }
2197             if (bORT || bODT)
2198             {
2199                 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2200                 fprintf(stderr, "End your selection with 0\n");
2201                 j     = -1;
2202                 orsel = NULL;
2203                 do
2204                 {
2205                     j++;
2206                     srenew(orsel, j+1);
2207                     if (1 != scanf("%d", &(orsel[j])))
2208                     {
2209                         gmx_fatal(FARGS, "Error reading user input");
2210                     }
2211                 }
2212                 while (orsel[j] > 0);
2213                 if (orsel[0] == -1)
2214                 {
2215                     fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2216                     norsel = nor;
2217                     srenew(orsel, nor);
2218                     for (i = 0; i < nor; i++)
2219                     {
2220                         orsel[i] = i;
2221                     }
2222                 }
2223                 else
2224                 {
2225                     /* Build the selection */
2226                     norsel = 0;
2227                     for (i = 0; i < j; i++)
2228                     {
2229                         for (k = 0; k < nor; k++)
2230                         {
2231                             if (or_label[k] == orsel[i])
2232                             {
2233                                 orsel[norsel] = k;
2234                                 norsel++;
2235                                 break;
2236                             }
2237                         }
2238                         if (k == nor)
2239                         {
2240                             fprintf(stderr, "Orientation restraint label %d not found\n",
2241                                     orsel[i]);
2242                         }
2243                     }
2244                 }
2245                 snew(odtleg, norsel);
2246                 for (i = 0; i < norsel; i++)
2247                 {
2248                     snew(odtleg[i], 256);
2249                     sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2250                 }
2251                 if (bORT)
2252                 {
2253                     fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2254                                     "Time (ps)", "", oenv);
2255                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2256                     {
2257                         fprintf(fort, "%s", orinst_sub);
2258                     }
2259                     xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2260                 }
2261                 if (bODT)
2262                 {
2263                     fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2264                                     "Orientation restraint deviation",
2265                                     "Time (ps)", "", oenv);
2266                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2267                     {
2268                         fprintf(fodt, "%s", orinst_sub);
2269                     }
2270                     xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2271                 }
2272             }
2273         }
2274         if (bOTEN)
2275         {
2276             foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2277                              "Order tensor", "Time (ps)", "", oenv);
2278             snew(otenleg, bOvec ? nex*12 : nex*3);
2279             for (i = 0; i < nex; i++)
2280             {
2281                 for (j = 0; j < 3; j++)
2282                 {
2283                     sprintf(buf, "eig%d", j+1);
2284                     otenleg[(bOvec ? 12 : 3)*i+j] = gmx_strdup(buf);
2285                 }
2286                 if (bOvec)
2287                 {
2288                     for (j = 0; j < 9; j++)
2289                     {
2290                         sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2291                         otenleg[12*i+3+j] = gmx_strdup(buf);
2292                     }
2293                 }
2294             }
2295             xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2296         }
2297     }
2298     else if (bDisRe)
2299     {
2300         nbounds = get_bounds(ftp2fn(efTPR, NFILE, fnm), &bounds, &index, &pair, &npairs,
2301                              &mtop, &top, &ir);
2302         snew(violaver, npairs);
2303         out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2304                        "Time (ps)", "nm", oenv);
2305         xvgr_legend(out, 2, drleg, oenv);
2306         if (bDRAll)
2307         {
2308             fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2309                                 "Time (ps)", "Distance (nm)", oenv);
2310             if (output_env_get_print_xvgr_codes(oenv))
2311             {
2312                 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2313                         ir.dr_tau);
2314             }
2315         }
2316     }
2317     else if (bDHDL)
2318     {
2319         get_dhdl_parms(ftp2fn(efTPR, NFILE, fnm), &ir);
2320     }
2321
2322     /* Initiate energies and set them to zero */
2323     edat.nsteps    = 0;
2324     edat.npoints   = 0;
2325     edat.nframes   = 0;
2326     edat.step      = NULL;
2327     edat.steps     = NULL;
2328     edat.points    = NULL;
2329     edat.bHaveSums = TRUE;
2330     snew(edat.s, nset);
2331
2332     /* Initiate counters */
2333     teller       = 0;
2334     teller_disre = 0;
2335     bFoundStart  = FALSE;
2336     start_step   = 0;
2337     start_t      = 0;
2338     do
2339     {
2340         /* This loop searches for the first frame (when -b option is given),
2341          * or when this has been found it reads just one energy frame
2342          */
2343         do
2344         {
2345             bCont = do_enx(fp, &(frame[NEXT]));
2346             if (bCont)
2347             {
2348                 timecheck = check_times(frame[NEXT].t);
2349             }
2350         }
2351         while (bCont && (timecheck < 0));
2352
2353         if ((timecheck == 0) && bCont)
2354         {
2355             /* We read a valid frame, so we can use it */
2356             fr = &(frame[NEXT]);
2357
2358             if (fr->nre > 0)
2359             {
2360                 /* The frame contains energies, so update cur */
2361                 cur  = NEXT;
2362
2363                 if (edat.nframes % 1000 == 0)
2364                 {
2365                     srenew(edat.step, edat.nframes+1000);
2366                     std::memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2367                     srenew(edat.steps, edat.nframes+1000);
2368                     std::memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2369                     srenew(edat.points, edat.nframes+1000);
2370                     std::memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2371
2372                     for (i = 0; i < nset; i++)
2373                     {
2374                         srenew(edat.s[i].ener, edat.nframes+1000);
2375                         std::memset(&(edat.s[i].ener[edat.nframes]), 0,
2376                                     1000*sizeof(edat.s[i].ener[0]));
2377                         srenew(edat.s[i].es, edat.nframes+1000);
2378                         std::memset(&(edat.s[i].es[edat.nframes]), 0,
2379                                     1000*sizeof(edat.s[i].es[0]));
2380                     }
2381                 }
2382
2383                 nfr            = edat.nframes;
2384                 edat.step[nfr] = fr->step;
2385
2386                 if (!bFoundStart)
2387                 {
2388                     bFoundStart = TRUE;
2389                     /* Initiate the previous step data */
2390                     start_step = fr->step;
2391                     start_t    = fr->t;
2392                     /* Initiate the energy sums */
2393                     edat.steps[nfr]  = 1;
2394                     edat.points[nfr] = 1;
2395                     for (i = 0; i < nset; i++)
2396                     {
2397                         sss                    = set[i];
2398                         edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2399                         edat.s[i].es[nfr].sum2 = 0;
2400                     }
2401                     edat.nsteps  = 1;
2402                     edat.npoints = 1;
2403                 }
2404                 else
2405                 {
2406                     edat.steps[nfr] = fr->nsteps;
2407
2408                     if (fr->nsum <= 1)
2409                     {
2410                         /* mdrun only calculated the energy at energy output
2411                          * steps. We don't need to check step intervals.
2412                          */
2413                         edat.points[nfr] = 1;
2414                         for (i = 0; i < nset; i++)
2415                         {
2416                             sss                    = set[i];
2417                             edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2418                             edat.s[i].es[nfr].sum2 = 0;
2419                         }
2420                         edat.npoints  += 1;
2421                         edat.bHaveSums = FALSE;
2422                     }
2423                     else if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2424                     {
2425                         /* We have statistics  to the previous frame */
2426                         edat.points[nfr] = fr->nsum;
2427                         for (i = 0; i < nset; i++)
2428                         {
2429                             sss                    = set[i];
2430                             edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2431                             edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2432                         }
2433                         edat.npoints += fr->nsum;
2434                     }
2435                     else
2436                     {
2437                         /* The interval does not match fr->nsteps:
2438                          * can not do exact averages.
2439                          */
2440                         edat.bHaveSums = FALSE;
2441                     }
2442
2443                     edat.nsteps = fr->step - start_step + 1;
2444                 }
2445                 for (i = 0; i < nset; i++)
2446                 {
2447                     edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2448                 }
2449             }
2450             /*
2451              * Define distance restraint legends. Can only be done after
2452              * the first frame has been read... (Then we know how many there are)
2453              */
2454             blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2455             if (bDisRe && bDRAll && !leg && blk_disre)
2456             {
2457                 t_iatom   *fa;
2458                 t_iparams *ip;
2459
2460                 fa = top->idef.il[F_DISRES].iatoms;
2461                 ip = top->idef.iparams;
2462                 if (blk_disre->nsub != 2 ||
2463                     (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2464                 {
2465                     gmx_incons("Number of disre sub-blocks not equal to 2");
2466                 }
2467
2468                 ndisre = blk_disre->sub[0].nr;
2469                 if (ndisre != top->idef.il[F_DISRES].nr/3)
2470                 {
2471                     gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2472                               ndisre, top->idef.il[F_DISRES].nr/3);
2473                 }
2474                 snew(pairleg, ndisre);
2475                 for (i = 0; i < ndisre; i++)
2476                 {
2477                     snew(pairleg[i], 30);
2478                     j = fa[3*i+1];
2479                     k = fa[3*i+2];
2480                     gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2481                     gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2482                     sprintf(pairleg[i], "%d %s %d %s (%d)",
2483                             resnr_j, anm_j, resnr_k, anm_k,
2484                             ip[fa[3*i]].disres.label);
2485                 }
2486                 set = select_it(ndisre, pairleg, &nset);
2487                 snew(leg, 2*nset);
2488                 for (i = 0; (i < nset); i++)
2489                 {
2490                     snew(leg[2*i], 32);
2491                     sprintf(leg[2*i],  "a %s", pairleg[set[i]]);
2492                     snew(leg[2*i+1], 32);
2493                     sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2494                 }
2495                 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2496             }
2497
2498             /*
2499              * Store energies for analysis afterwards...
2500              */
2501             if (!bDisRe && !bDHDL && (fr->nre > 0))
2502             {
2503                 if (edat.nframes % 1000 == 0)
2504                 {
2505                     srenew(time, edat.nframes+1000);
2506                 }
2507                 time[edat.nframes] = fr->t;
2508                 edat.nframes++;
2509             }
2510             /*
2511              * Printing time, only when we do not want to skip frames
2512              */
2513             if (!skip || teller % skip == 0)
2514             {
2515                 if (bDisRe)
2516                 {
2517                     /*******************************************
2518                      * D I S T A N C E   R E S T R A I N T S
2519                      *******************************************/
2520                     if (ndisre > 0)
2521                     {
2522                         GMX_RELEASE_ASSERT(blk_disre != NULL, "Trying to dereference NULL blk_disre pointer");
2523  #ifndef GMX_DOUBLE
2524                         float  *disre_rt     =     blk_disre->sub[0].fval;
2525                         float  *disre_rm3tav = blk_disre->sub[1].fval;
2526  #else
2527                         double *disre_rt     =     blk_disre->sub[0].dval;
2528                         double *disre_rm3tav = blk_disre->sub[1].dval;
2529  #endif
2530
2531                         print_time(out, fr->t);
2532                         if (violaver == NULL)
2533                         {
2534                             snew(violaver, ndisre);
2535                         }
2536
2537                         /* Subtract bounds from distances, to calculate violations */
2538                         calc_violations(disre_rt, disre_rm3tav,
2539                                         nbounds, pair, bounds, violaver, &sumt, &sumaver);
2540
2541                         fprintf(out, "  %8.4f  %8.4f\n", sumaver, sumt);
2542                         if (bDRAll)
2543                         {
2544                             print_time(fp_pairs, fr->t);
2545                             for (i = 0; (i < nset); i++)
2546                             {
2547                                 sss = set[i];
2548                                 fprintf(fp_pairs, "  %8.4f", mypow(disre_rm3tav[sss], minthird));
2549                                 fprintf(fp_pairs, "  %8.4f", disre_rt[sss]);
2550                             }
2551                             fprintf(fp_pairs, "\n");
2552                         }
2553                         teller_disre++;
2554                     }
2555                 }
2556                 else if (bDHDL)
2557                 {
2558                     do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2559                 }
2560
2561                 /*******************************************
2562                  * E N E R G I E S
2563                  *******************************************/
2564                 else
2565                 {
2566                     if (fr->nre > 0)
2567                     {
2568                         if (bPrAll)
2569                         {
2570                             /* We skip frames with single points (usually only the first frame),
2571                              * since they would result in an average plot with outliers.
2572                              */
2573                             if (fr->nsum > 1)
2574                             {
2575                                 print_time(out, fr->t);
2576                                 print1(out, bDp, fr->ener[set[0]].e);
2577                                 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2578                                 print1(out, bDp, std::sqrt(fr->ener[set[0]].eav/fr->nsum));
2579                                 fprintf(out, "\n");
2580                             }
2581                         }
2582                         else
2583                         {
2584                             print_time(out, fr->t);
2585                             if (bSum)
2586                             {
2587                                 sum = 0;
2588                                 for (i = 0; i < nset; i++)
2589                                 {
2590                                     sum += fr->ener[set[i]].e;
2591                                 }
2592                                 print1(out, bDp, sum/nmol-ezero);
2593                             }
2594                             else
2595                             {
2596                                 for (i = 0; (i < nset); i++)
2597                                 {
2598                                     if (bIsEner[i])
2599                                     {
2600                                         print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2601                                     }
2602                                     else
2603                                     {
2604                                         print1(out, bDp, fr->ener[set[i]].e);
2605                                     }
2606                                 }
2607                             }
2608                             fprintf(out, "\n");
2609                         }
2610                     }
2611                     blk = find_block_id_enxframe(fr, enx_i, NULL);
2612                     if (bORIRE && blk)
2613                     {
2614 #ifndef GMX_DOUBLE
2615                         xdr_datatype dt = xdr_datatype_float;
2616 #else
2617                         xdr_datatype dt = xdr_datatype_double;
2618 #endif
2619                         real        *vals;
2620
2621                         if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2622                         {
2623                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2624                         }
2625 #ifndef GMX_DOUBLE
2626                         vals = blk->sub[0].fval;
2627 #else
2628                         vals = blk->sub[0].dval;
2629 #endif
2630
2631                         if (blk->sub[0].nr != nor)
2632                         {
2633                             gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2634                         }
2635                         if (bORA || bODA)
2636                         {
2637                             for (i = 0; i < nor; i++)
2638                             {
2639                                 orient[i] += vals[i];
2640                             }
2641                         }
2642                         if (bODR)
2643                         {
2644                             for (i = 0; i < nor; i++)
2645                             {
2646                                 odrms[i] += sqr(vals[i]-oobs[i]);
2647                             }
2648                         }
2649                         if (bORT)
2650                         {
2651                             fprintf(fort, "  %10f", fr->t);
2652                             for (i = 0; i < norsel; i++)
2653                             {
2654                                 fprintf(fort, " %g", vals[orsel[i]]);
2655                             }
2656                             fprintf(fort, "\n");
2657                         }
2658                         if (bODT)
2659                         {
2660                             fprintf(fodt, "  %10f", fr->t);
2661                             for (i = 0; i < norsel; i++)
2662                             {
2663                                 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2664                             }
2665                             fprintf(fodt, "\n");
2666                         }
2667                         norfr++;
2668                     }
2669                     blk = find_block_id_enxframe(fr, enxORT, NULL);
2670                     if (bOTEN && blk)
2671                     {
2672 #ifndef GMX_DOUBLE
2673                         xdr_datatype dt = xdr_datatype_float;
2674 #else
2675                         xdr_datatype dt = xdr_datatype_double;
2676 #endif
2677                         real        *vals;
2678
2679                         if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2680                         {
2681                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2682                         }
2683 #ifndef GMX_DOUBLE
2684                         vals = blk->sub[0].fval;
2685 #else
2686                         vals = blk->sub[0].dval;
2687 #endif
2688
2689                         if (blk->sub[0].nr != nex*12)
2690                         {
2691                             gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2692                                       blk->sub[0].nr/12, nex);
2693                         }
2694                         fprintf(foten, "  %10f", fr->t);
2695                         for (i = 0; i < nex; i++)
2696                         {
2697                             for (j = 0; j < (bOvec ? 12 : 3); j++)
2698                             {
2699                                 fprintf(foten, " %g", vals[i*12+j]);
2700                             }
2701                         }
2702                         fprintf(foten, "\n");
2703                     }
2704                 }
2705             }
2706             teller++;
2707         }
2708     }
2709     while (bCont && (timecheck == 0));
2710
2711     fprintf(stderr, "\n");
2712     close_enx(fp);
2713     if (out)
2714     {
2715         xvgrclose(out);
2716     }
2717
2718     if (bDRAll)
2719     {
2720         xvgrclose(fp_pairs);
2721     }
2722
2723     if (bORT)
2724     {
2725         xvgrclose(fort);
2726     }
2727     if (bODT)
2728     {
2729         xvgrclose(fodt);
2730     }
2731     if (bORA)
2732     {
2733         out = xvgropen(opt2fn("-ora", NFILE, fnm),
2734                        "Average calculated orientations",
2735                        "Restraint label", "", oenv);
2736         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2737         {
2738             fprintf(out, "%s", orinst_sub);
2739         }
2740         for (i = 0; i < nor; i++)
2741         {
2742             fprintf(out, "%5d  %g\n", or_label[i], orient[i]/norfr);
2743         }
2744         xvgrclose(out);
2745     }
2746     if (bODA)
2747     {
2748         out = xvgropen(opt2fn("-oda", NFILE, fnm),
2749                        "Average restraint deviation",
2750                        "Restraint label", "", oenv);
2751         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2752         {
2753             fprintf(out, "%s", orinst_sub);
2754         }
2755         for (i = 0; i < nor; i++)
2756         {
2757             fprintf(out, "%5d  %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2758         }
2759         xvgrclose(out);
2760     }
2761     if (bODR)
2762     {
2763         out = xvgropen(opt2fn("-odr", NFILE, fnm),
2764                        "RMS orientation restraint deviations",
2765                        "Restraint label", "", oenv);
2766         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2767         {
2768             fprintf(out, "%s", orinst_sub);
2769         }
2770         for (i = 0; i < nor; i++)
2771         {
2772             fprintf(out, "%5d  %g\n", or_label[i], std::sqrt(odrms[i]/norfr));
2773         }
2774         xvgrclose(out);
2775     }
2776     if (bOTEN)
2777     {
2778         xvgrclose(foten);
2779     }
2780
2781     if (bDisRe)
2782     {
2783         analyse_disre(opt2fn("-viol", NFILE, fnm),
2784                       teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2785     }
2786     else if (bDHDL)
2787     {
2788         if (fp_dhdl)
2789         {
2790             gmx_fio_fclose(fp_dhdl);
2791             printf("\n\nWrote %d lambda values with %d samples as ",
2792                    dh_lambdas, dh_samples);
2793             if (dh_hists > 0)
2794             {
2795                 printf("%d dH histograms ", dh_hists);
2796             }
2797             if (dh_blocks > 0)
2798             {
2799                 printf("%d dH data blocks ", dh_blocks);
2800             }
2801             printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2802
2803         }
2804         else
2805         {
2806             gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2807         }
2808
2809     }
2810     else
2811     {
2812         double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2813         analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2814                      bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa),
2815                      bVisco, opt2fn("-vis", NFILE, fnm),
2816                      nmol,
2817                      start_step, start_t, frame[cur].step, frame[cur].t,
2818                      reftemp, &edat,
2819                      nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2820                      oenv);
2821         if (bFluctProps)
2822         {
2823             calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2824                                    nbmin, nbmax);
2825         }
2826     }
2827     if (opt2bSet("-f2", NFILE, fnm))
2828     {
2829         fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2830             reftemp, nset, set, leg, &edat, time, oenv);
2831     }
2832
2833     {
2834         const char *nxy = "-nxy";
2835
2836         do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2837         do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2838         do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2839         do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2840         do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2841         do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2842         do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2843         do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2844         do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);
2845     }
2846
2847     return 0;
2848 }