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