src/gromacs/gmxana/gmx_energy.c

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