src/gromacs/gmxana/gmx_energy.c

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