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