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 "gromacs/fileio/enxio.h"
  50 #include "gromacs/commandline/pargs.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 "gromacs/fileio/tpxio.h"
  58 #include "gromacs/fileio/trxio.h"
  59 #include "viewit.h"
  60 #include "mtop_util.h"
  61 #include "gmx_ana.h"
  62 #include "mdebin.h"
  63
  64 static real       minthird = -1.0/3.0, minsixth = -1.0/6.0;
  65
  66 typedef struct {
  67     real sum;
  68     real sum2;
  69 } exactsum_t;
  70
  71 typedef struct {
  72     real       *ener;
  73     exactsum_t *es;
  74     gmx_bool    bExactStat;
  75     double      av;
  76     double      rmsd;
  77     double      ee;
  78     double      slope;
  79 } enerdat_t;
  80
  81 typedef struct {
  82     gmx_int64_t      nsteps;
  83     gmx_int64_t      npoints;
  84     int              nframes;
  85     int             *step;
  86     int             *steps;
  87     int             *points;
  88     enerdat_t       *s;
  89 } enerdata_t;
  90
  91 static double mypow(double x, double y)
  92 {
  93     if (x > 0)
  94     {
  95         return pow(x, y);
  96     }
  97     else
  98     {
  99         return 0.0;
 100     }
 101 }
 102
 103 static int *select_it(int nre, char *nm[], int *nset)
 104 {
 105     gmx_bool *bE;
 106     int       n, k, j, i;
 107     int      *set;
 108     gmx_bool  bVerbose = TRUE;
 109
 110     if ((getenv("VERBOSE")) != NULL)
 111     {
 112         bVerbose = FALSE;
 113     }
 114
 115     fprintf(stderr, "Select the terms you want from the following list\n");
 116     fprintf(stderr, "End your selection with 0\n");
 117
 118     if (bVerbose)
 119     {
 120         for (k = 0; (k < nre); )
 121         {
 122             for (j = 0; (j < 4) && (k < nre); j++, k++)
 123             {
 124                 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
 125             }
 126             fprintf(stderr, "\n");
 127         }
 128     }
 129
 130     snew(bE, nre);
 131     do
 132     {
 133         if (1 != scanf("%d", &n))
 134         {
 135             gmx_fatal(FARGS, "Error reading user input");
 136         }
 137         if ((n > 0) && (n <= nre))
 138         {
 139             bE[n-1] = TRUE;
 140         }
 141     }
 142     while (n != 0);
 143
 144     snew(set, nre);
 145     for (i = (*nset) = 0; (i < nre); i++)
 146     {
 147         if (bE[i])
 148         {
 149             set[(*nset)++] = i;
 150         }
 151     }
 152
 153     sfree(bE);
 154
 155     return set;
 156 }
 157
 158 static void chomp(char *buf)
 159 {
 160     int len = strlen(buf);
 161
 162     while ((len > 0) && (buf[len-1] == '\n'))
 163     {
 164         buf[len-1] = '\0';
 165         len--;
 166     }
 167 }
 168
 169 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
 170 {
 171     gmx_bool   *bE;
 172     int         n, k, kk, j, i, nmatch, nind, nss;
 173     int        *set;
 174     gmx_bool    bEOF, bVerbose = TRUE, bLong = FALSE;
 175     char       *ptr, buf[STRLEN];
 176     const char *fm4   = "%3d  %-14s";
 177     const char *fm2   = "%3d  %-34s";
 178     char      **newnm = NULL;
 179
 180     if ((getenv("VERBOSE")) != NULL)
 181     {
 182         bVerbose = FALSE;
 183     }
 184
 185     fprintf(stderr, "\n");
 186     fprintf(stderr, "Select the terms you want from the following list by\n");
 187     fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
 188     fprintf(stderr, "End your selection with an empty line or a zero.\n");
 189     fprintf(stderr, "-------------------------------------------------------------------\n");
 190
 191     snew(newnm, nre);
 192     j = 0;
 193     for (k = 0; k < nre; k++)
 194     {
 195         newnm[k] = strdup(nm[k].name);
 196         /* Insert dashes in all the names */
 197         while ((ptr = strchr(newnm[k], ' ')) != NULL)
 198         {
 199             *ptr = '-';
 200         }
 201         if (bVerbose)
 202         {
 203             if (j == 0)
 204             {
 205                 if (k > 0)
 206                 {
 207                     fprintf(stderr, "\n");
 208                 }
 209                 bLong = FALSE;
 210                 for (kk = k; kk < k+4; kk++)
 211                 {
 212                     if (kk < nre && strlen(nm[kk].name) > 14)
 213                     {
 214                         bLong = TRUE;
 215                     }
 216                 }
 217             }
 218             else
 219             {
 220                 fprintf(stderr, " ");
 221             }
 222             if (!bLong)
 223             {
 224                 fprintf(stderr, fm4, k+1, newnm[k]);
 225                 j++;
 226                 if (j == 4)
 227                 {
 228                     j = 0;
 229                 }
 230             }
 231             else
 232             {
 233                 fprintf(stderr, fm2, k+1, newnm[k]);
 234                 j++;
 235                 if (j == 2)
 236                 {
 237                     j = 0;
 238                 }
 239             }
 240         }
 241     }
 242     if (bVerbose)
 243     {
 244         fprintf(stderr, "\n\n");
 245     }
 246
 247     snew(bE, nre);
 248
 249     bEOF = FALSE;
 250     while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
 251     {
 252         /* Remove newlines */
 253         chomp(buf);
 254
 255         /* Remove spaces */
 256         trim(buf);
 257
 258         /* Empty line means end of input */
 259         bEOF = (strlen(buf) == 0);
 260         if (!bEOF)
 261         {
 262             ptr = buf;
 263             do
 264             {
 265                 if (!bEOF)
 266                 {
 267                     /* First try to read an integer */
 268                     nss   = sscanf(ptr, "%d", &nind);
 269                     if (nss == 1)
 270                     {
 271                         /* Zero means end of input */
 272                         if (nind == 0)
 273                         {
 274                             bEOF = TRUE;
 275                         }
 276                         else if ((1 <= nind) && (nind <= nre))
 277                         {
 278                             bE[nind-1] = TRUE;
 279                         }
 280                         else
 281                         {
 282                             fprintf(stderr, "number %d is out of range\n", nind);
 283                         }
 284                     }
 285                     else
 286                     {
 287                         /* Now try to read a string */
 288                         i      = strlen(ptr);
 289                         nmatch = 0;
 290                         for (nind = 0; nind < nre; nind++)
 291                         {
 292                             if (gmx_strcasecmp(newnm[nind], ptr) == 0)
 293                             {
 294                                 bE[nind] = TRUE;
 295                                 nmatch++;
 296                             }
 297                         }
 298                         if (nmatch == 0)
 299                         {
 300                             i      = strlen(ptr);
 301                             nmatch = 0;
 302                             for (nind = 0; nind < nre; nind++)
 303                             {
 304                                 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
 305                                 {
 306                                     bE[nind] = TRUE;
 307                                     nmatch++;
 308                                 }
 309                             }
 310                             if (nmatch == 0)
 311                             {
 312                                 fprintf(stderr, "String '%s' does not match anything\n", ptr);
 313                             }
 314                         }
 315                     }
 316                 }
 317                 /* Look for the first space, and remove spaces from there */
 318                 if ((ptr = strchr(ptr, ' ')) != NULL)
 319                 {
 320                     trim(ptr);
 321                 }
 322             }
 323             while (!bEOF && (ptr && (strlen(ptr) > 0)));
 324         }
 325     }
 326
 327     snew(set, nre);
 328     for (i = (*nset) = 0; (i < nre); i++)
 329     {
 330         if (bE[i])
 331         {
 332             set[(*nset)++] = i;
 333         }
 334     }
 335
 336     sfree(bE);
 337
 338     if (*nset == 0)
 339     {
 340         gmx_fatal(FARGS, "No energy terms selected");
 341     }
 342
 343     for (i = 0; (i < nre); i++)
 344     {
 345         sfree(newnm[i]);
 346     }
 347     sfree(newnm);
 348
 349     return set;
 350 }
 351
 352 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
 353 {
 354     gmx_mtop_t  mtop;
 355     int         natoms;
 356     t_iatom    *iatom;
 357     matrix      box;
 358
 359     /* all we need is the ir to be able to write the label */
 360     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
 361 }
 362
 363 static void get_orires_parms(const char *topnm,
 364                              int *nor, int *nex, int **label, real **obs)
 365 {
 366     gmx_mtop_t      mtop;
 367     gmx_localtop_t *top;
 368     t_inputrec      ir;
 369     t_iparams      *ip;
 370     int             natoms, i;
 371     t_iatom        *iatom;
 372     int             nb;
 373     matrix          box;
 374
 375     read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
 376     top = gmx_mtop_generate_local_top(&mtop, &ir);
 377
 378     ip       = top->idef.iparams;
 379     iatom    = top->idef.il[F_ORIRES].iatoms;
 380
 381     /* Count how many distance restraint there are... */
 382     nb = top->idef.il[F_ORIRES].nr;
 383     if (nb == 0)
 384     {
 385         gmx_fatal(FARGS, "No orientation restraints in topology!\n");
 386     }
 387
 388     *nor = nb/3;
 389     *nex = 0;
 390     snew(*label, *nor);
 391     snew(*obs, *nor);
 392     for (i = 0; i < nb; i += 3)
 393     {
 394         (*label)[i/3] = ip[iatom[i]].orires.label;
 395         (*obs)[i/3]   = ip[iatom[i]].orires.obs;
 396         if (ip[iatom[i]].orires.ex >= *nex)
 397         {
 398             *nex = ip[iatom[i]].orires.ex+1;
 399         }
 400     }
 401     fprintf(stderr, "Found %d orientation restraints with %d experiments",
 402             *nor, *nex);
 403 }
 404
 405 static int get_bounds(const char *topnm,
 406                       real **bounds, int **index, int **dr_pair, int *npairs,
 407                       gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
 408 {
 409     gmx_localtop_t *top;
 410     t_functype     *functype;
 411     t_iparams      *ip;
 412     int             natoms, i, j, k, type, ftype, natom;
 413     t_ilist        *disres;
 414     t_iatom        *iatom;
 415     real           *b;
 416     int            *ind, *pair;
 417     int             nb, label1;
 418     matrix          box;
 419
 420     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
 421     snew(*ltop, 1);
 422     top   = gmx_mtop_generate_local_top(mtop, ir);
 423     *ltop = top;
 424
 425     functype = top->idef.functype;
 426     ip       = top->idef.iparams;
 427
 428     /* Count how many distance restraint there are... */
 429     nb = top->idef.il[F_DISRES].nr;
 430     if (nb == 0)
 431     {
 432         gmx_fatal(FARGS, "No distance restraints in topology!\n");
 433     }
 434
 435     /* Allocate memory */
 436     snew(b, nb);
 437     snew(ind, nb);
 438     snew(pair, nb+1);
 439
 440     /* Fill the bound array */
 441     nb = 0;
 442     for (i = 0; (i < top->idef.ntypes); i++)
 443     {
 444         ftype = functype[i];
 445         if (ftype == F_DISRES)
 446         {
 447
 448             label1  = ip[i].disres.label;
 449             b[nb]   = ip[i].disres.up1;
 450             ind[nb] = label1;
 451             nb++;
 452         }
 453     }
 454     *bounds = b;
 455
 456     /* Fill the index array */
 457     label1  = -1;
 458     disres  = &(top->idef.il[F_DISRES]);
 459     iatom   = disres->iatoms;
 460     for (i = j = k = 0; (i < disres->nr); )
 461     {
 462         type  = iatom[i];
 463         ftype = top->idef.functype[type];
 464         natom = interaction_function[ftype].nratoms+1;
 465         if (label1 != top->idef.iparams[type].disres.label)
 466         {
 467             pair[j] = k;
 468             label1  = top->idef.iparams[type].disres.label;
 469             j++;
 470         }
 471         k++;
 472         i += natom;
 473     }
 474     pair[j]  = k;
 475     *npairs  = k;
 476     if (j != nb)
 477     {
 478         gmx_incons("get_bounds for distance restraints");
 479     }
 480
 481     *index   = ind;
 482     *dr_pair = pair;
 483
 484     return nb;
 485 }
 486
 487 static void calc_violations(real rt[], real rav3[], int nb, int index[],
 488                             real bounds[], real *viol, double *st, double *sa)
 489 {
 490     const   real sixth = 1.0/6.0;
 491     int          i, j;
 492     double       rsum, rav, sumaver, sumt;
 493
 494     sumaver = 0;
 495     sumt    = 0;
 496     for (i = 0; (i < nb); i++)
 497     {
 498         rsum = 0.0;
 499         rav  = 0.0;
 500         for (j = index[i]; (j < index[i+1]); j++)
 501         {
 502             if (viol)
 503             {
 504                 viol[j] += mypow(rt[j], -3.0);
 505             }
 506             rav     += sqr(rav3[j]);
 507             rsum    += mypow(rt[j], -6);
 508         }
 509         rsum    = max(0.0, mypow(rsum, -sixth)-bounds[i]);
 510         rav     = max(0.0, mypow(rav, -sixth)-bounds[i]);
 511
 512         sumt    += rsum;
 513         sumaver += rav;
 514     }
 515     *st = sumt;
 516     *sa = sumaver;
 517 }
 518
 519 static void analyse_disre(const char *voutfn,    int nframes,
 520                           real violaver[], real bounds[], int index[],
 521                           int pair[],      int nbounds,
 522                           const output_env_t oenv)
 523 {
 524     FILE   *vout;
 525     double  sum, sumt, sumaver;
 526     int     i, j;
 527
 528     /* Subtract bounds from distances, to calculate violations */
 529     calc_violations(violaver, violaver,
 530                     nbounds, pair, bounds, NULL, &sumt, &sumaver);
 531
 532 #ifdef DEBUG
 533     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 534             sumaver);
 535     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
 536             sumt);
 537 #endif
 538     vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
 539                     oenv);
 540     sum  = 0.0;
 541     sumt = 0.0;
 542     for (i = 0; (i < nbounds); i++)
 543     {
 544         /* Do ensemble averaging */
 545         sumaver = 0;
 546         for (j = pair[i]; (j < pair[i+1]); j++)
 547         {
 548             sumaver += sqr(violaver[j]/nframes);
 549         }
 550         sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
 551
 552         sumt   += sumaver;
 553         sum     = max(sum, sumaver);
 554         fprintf(vout, "%10d  %10.5e\n", index[i], sumaver);
 555     }
 556 #ifdef DEBUG
 557     for (j = 0; (j < dr.ndr); j++)
 558     {
 559         fprintf(vout, "%10d  %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
 560     }
 561 #endif
 562     ffclose(vout);
 563
 564     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 565             sumt);
 566     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
 567
 568     do_view(oenv, voutfn, "-graphtype bar");
 569 }
 570
 571 static void einstein_visco(const char *fn, const char *fni, int nsets,
 572                            int nframes, real **sum,
 573                            real V, real T, int nsteps, double time[],
 574                            const output_env_t oenv)
 575 {
 576     FILE  *fp0, *fp1;
 577     double av[4], avold[4];
 578     double fac, dt, di;
 579     int    i, j, m, nf4;
 580
 581     if (nframes < 1)
 582     {
 583         return;
 584     }
 585
 586     dt  = (time[1]-time[0]);
 587     nf4 = nframes/4+1;
 588
 589     for (i = 0; i <= nsets; i++)
 590     {
 591         avold[i] = 0;
 592     }
 593     fp0 = xvgropen(fni, "Shear viscosity integral",
 594                    "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
 595     fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
 596                    "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
 597     for (i = 1; i < nf4; i++)
 598     {
 599         fac = dt*nframes/nsteps;
 600         for (m = 0; m <= nsets; m++)
 601         {
 602             av[m] = 0;
 603         }
 604         for (j = 0; j < nframes-i; j++)
 605         {
 606             for (m = 0; m < nsets; m++)
 607             {
 608                 di   = sqr(fac*(sum[m][j+i]-sum[m][j]));
 609
 610                 av[m]     += di;
 611                 av[nsets] += di/nsets;
 612             }
 613         }
 614         /* Convert to SI for the viscosity */
 615         fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
 616         fprintf(fp0, "%10g", time[i]-time[0]);
 617         for (m = 0; (m <= nsets); m++)
 618         {
 619             av[m] = fac*av[m];
 620             fprintf(fp0, "  %10g", av[m]);
 621         }
 622         fprintf(fp0, "\n");
 623         fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]);
 624         for (m = 0; (m <= nsets); m++)
 625         {
 626             fprintf(fp1, "  %10g", (av[m]-avold[m])/dt);
 627             avold[m] = av[m];
 628         }
 629         fprintf(fp1, "\n");
 630     }
 631     ffclose(fp0);
 632     ffclose(fp1);
 633 }
 634
 635 typedef struct {
 636     gmx_int64_t     np;
 637     double          sum;
 638     double          sav;
 639     double          sav2;
 640 } ee_sum_t;
 641
 642 typedef struct {
 643     int             b;
 644     ee_sum_t        sum;
 645     gmx_int64_t     nst;
 646     gmx_int64_t     nst_min;
 647 } ener_ee_t;
 648
 649 static void clear_ee_sum(ee_sum_t *ees)
 650 {
 651     ees->sav  = 0;
 652     ees->sav2 = 0;
 653     ees->np   = 0;
 654     ees->sum  = 0;
 655 }
 656
 657 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
 658 {
 659     ees->np  += np;
 660     ees->sum += sum;
 661 }
 662
 663 static void add_ee_av(ee_sum_t *ees)
 664 {
 665     double av;
 666
 667     av         = ees->sum/ees->np;
 668     ees->sav  += av;
 669     ees->sav2 += av*av;
 670     ees->np    = 0;
 671     ees->sum   = 0;
 672 }
 673
 674 static double calc_ee2(int nb, ee_sum_t *ees)
 675 {
 676     return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
 677 }
 678
 679 static void set_ee_av(ener_ee_t *eee)
 680 {
 681     if (debug)
 682     {
 683         char buf[STEPSTRSIZE];
 684         fprintf(debug, "Storing average for err.est.: %s steps\n",
 685                 gmx_step_str(eee->nst, buf));
 686     }
 687     add_ee_av(&eee->sum);
 688     eee->b++;
 689     if (eee->b == 1 || eee->nst < eee->nst_min)
 690     {
 691         eee->nst_min = eee->nst;
 692     }
 693     eee->nst = 0;
 694 }
 695
 696 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
 697 {
 698     int             nb, i, f, nee;
 699     double          sum, sum2, sump, see2;
 700     gmx_int64_t     steps, np, p, bound_nb;
 701     enerdat_t      *ed;
 702     exactsum_t     *es;
 703     gmx_bool        bAllZero;
 704     double          x, sx, sy, sxx, sxy;
 705     ener_ee_t      *eee;
 706
 707     /* Check if we have exact statistics over all points */
 708     for (i = 0; i < nset; i++)
 709     {
 710         ed             = &edat->s[i];
 711         ed->bExactStat = FALSE;
 712         if (edat->npoints > 0)
 713         {
 714             /* All energy file sum entries 0 signals no exact sums.
 715              * But if all energy values are 0, we still have exact sums.
 716              */
 717             bAllZero = TRUE;
 718             for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
 719             {
 720                 if (ed->ener[i] != 0)
 721                 {
 722                     bAllZero = FALSE;
 723                 }
 724                 ed->bExactStat = (ed->es[f].sum != 0);
 725             }
 726             if (bAllZero)
 727             {
 728                 ed->bExactStat = TRUE;
 729             }
 730         }
 731     }
 732
 733     snew(eee, nbmax+1);
 734     for (i = 0; i < nset; i++)
 735     {
 736         ed = &edat->s[i];
 737
 738         sum  = 0;
 739         sum2 = 0;
 740         np   = 0;
 741         sx   = 0;
 742         sy   = 0;
 743         sxx  = 0;
 744         sxy  = 0;
 745         for (nb = nbmin; nb <= nbmax; nb++)
 746         {
 747             eee[nb].b     = 0;
 748             clear_ee_sum(&eee[nb].sum);
 749             eee[nb].nst     = 0;
 750             eee[nb].nst_min = 0;
 751         }
 752         for (f = 0; f < edat->nframes; f++)
 753         {
 754             es = &ed->es[f];
 755
 756             if (ed->bExactStat)
 757             {
 758                 /* Add the sum and the sum of variances to the totals. */
 759                 p     = edat->points[f];
 760                 sump  = es->sum;
 761                 sum2 += es->sum2;
 762                 if (np > 0)
 763                 {
 764                     sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
 765                         *np*(np + p)/p;
 766                 }
 767             }
 768             else
 769             {
 770                 /* Add a single value to the sum and sum of squares. */
 771                 p     = 1;
 772                 sump  = ed->ener[f];
 773                 sum2 += dsqr(sump);
 774             }
 775
 776             /* sum has to be increased after sum2 */
 777             np  += p;
 778             sum += sump;
 779
 780             /* For the linear regression use variance 1/p.
 781              * Note that sump is the sum, not the average, so we don't need p*.
 782              */
 783             x    = edat->step[f] - 0.5*(edat->steps[f] - 1);
 784             sx  += p*x;
 785             sy  += sump;
 786             sxx += p*x*x;
 787             sxy += x*sump;
 788
 789             for (nb = nbmin; nb <= nbmax; nb++)
 790             {
 791                 /* Check if the current end step is closer to the desired
 792                  * block boundary than the next end step.
 793                  */
 794                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 795                 if (eee[nb].nst > 0 &&
 796                     bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
 797                 {
 798                     set_ee_av(&eee[nb]);
 799                 }
 800                 if (f == 0)
 801                 {
 802                     eee[nb].nst = 1;
 803                 }
 804                 else
 805                 {
 806                     eee[nb].nst += edat->step[f] - edat->step[f-1];
 807                 }
 808                 if (ed->bExactStat)
 809                 {
 810                     add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
 811                 }
 812                 else
 813                 {
 814                     add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
 815                 }
 816                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 817                 if (edat->step[f]*nb >= bound_nb)
 818                 {
 819                     set_ee_av(&eee[nb]);
 820                 }
 821             }
 822         }
 823
 824         edat->s[i].av = sum/np;
 825         if (ed->bExactStat)
 826         {
 827             edat->s[i].rmsd = sqrt(sum2/np);
 828         }
 829         else
 830         {
 831             edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
 832         }
 833
 834         if (edat->nframes > 1)
 835         {
 836             edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
 837         }
 838         else
 839         {
 840             edat->s[i].slope = 0;
 841         }
 842
 843         nee  = 0;
 844         see2 = 0;
 845         for (nb = nbmin; nb <= nbmax; nb++)
 846         {
 847             /* Check if we actually got nb blocks and if the smallest
 848              * block is not shorter than 80% of the average.
 849              */
 850             if (debug)
 851             {
 852                 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
 853                 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
 854                         nb, eee[nb].b,
 855                         gmx_step_str(eee[nb].nst_min, buf1),
 856                         gmx_step_str(edat->nsteps, buf2));
 857             }
 858             if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
 859             {
 860                 see2 += calc_ee2(nb, &eee[nb].sum);
 861                 nee++;
 862             }
 863         }
 864         if (nee > 0)
 865         {
 866             edat->s[i].ee = sqrt(see2/nee);
 867         }
 868         else
 869         {
 870             edat->s[i].ee = -1;
 871         }
 872     }
 873     sfree(eee);
 874 }
 875
 876 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
 877 {
 878     enerdata_t *esum;
 879     enerdat_t  *s;
 880     int         f, i;
 881     double      sum;
 882
 883     snew(esum, 1);
 884     *esum = *edat;
 885     snew(esum->s, 1);
 886     s = &esum->s[0];
 887     snew(s->ener, esum->nframes);
 888     snew(s->es, esum->nframes);
 889
 890     s->bExactStat = TRUE;
 891     s->slope      = 0;
 892     for (i = 0; i < nset; i++)
 893     {
 894         if (!edat->s[i].bExactStat)
 895         {
 896             s->bExactStat = FALSE;
 897         }
 898         s->slope += edat->s[i].slope;
 899     }
 900
 901     for (f = 0; f < edat->nframes; f++)
 902     {
 903         sum = 0;
 904         for (i = 0; i < nset; i++)
 905         {
 906             sum += edat->s[i].ener[f];
 907         }
 908         s->ener[f] = sum;
 909         sum        = 0;
 910         for (i = 0; i < nset; i++)
 911         {
 912             sum += edat->s[i].es[f].sum;
 913         }
 914         s->es[f].sum  = sum;
 915         s->es[f].sum2 = 0;
 916     }
 917
 918     calc_averages(1, esum, nbmin, nbmax);
 919
 920     return esum;
 921 }
 922
 923 static char *ee_pr(double ee, char *buf)
 924 {
 925     char   tmp[100];
 926     double rnd;
 927
 928     if (ee < 0)
 929     {
 930         sprintf(buf, "%s", "--");
 931     }
 932     else
 933     {
 934         /* Round to two decimals by printing. */
 935         sprintf(tmp, "%.1e", ee);
 936         sscanf(tmp, "%lf", &rnd);
 937         sprintf(buf, "%g", rnd);
 938     }
 939
 940     return buf;
 941 }
 942
 943 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
 944 {
 945 /* Remove the drift by performing a fit to y = ax+b.
 946    Uses 5 iterations. */
 947     int    i, j, k;
 948     double delta, d, sd, sd2;
 949
 950     edat->npoints = edat->nframes;
 951     edat->nsteps  = edat->nframes;
 952
 953     for (k = 0; (k < 5); k++)
 954     {
 955         for (i = 0; (i < nset); i++)
 956         {
 957             delta = edat->s[i].slope*dt;
 958
 959             if (NULL != debug)
 960             {
 961                 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
 962             }
 963
 964             for (j = 0; (j < edat->nframes); j++)
 965             {
 966                 edat->s[i].ener[j]   -= j*delta;
 967                 edat->s[i].es[j].sum  = 0;
 968                 edat->s[i].es[j].sum2 = 0;
 969             }
 970         }
 971         calc_averages(nset, edat, nbmin, nbmax);
 972     }
 973 }
 974
 975 static void calc_fluctuation_props(FILE *fp,
 976                                    gmx_bool bDriftCorr, real dt,
 977                                    int nset, int nmol,
 978                                    char **leg, enerdata_t *edat,
 979                                    int nbmin, int nbmax)
 980 {
 981     int    i, j;
 982     double vvhh, vv, v, h, hh2, vv2, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
 983     double NANO3;
 984     enum {
 985         eVol, eEnth, eTemp, eEtot, eNR
 986     };
 987     const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
 988     int         ii[eNR];
 989
 990     NANO3 = NANO*NANO*NANO;
 991     if (!bDriftCorr)
 992     {
 993         fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
 994                 "for spurious drift in the graphs. Note that this is not\n"
 995                 "a substitute for proper equilibration and sampling!\n");
 996     }
 997     else
 998     {
 999         remove_drift(nset, nbmin, nbmax, dt, edat);
1000     }
1001     for (i = 0; (i < eNR); i++)
1002     {
1003         for (ii[i] = 0; (ii[i] < nset &&
1004                          (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
1005         {
1006             ;
1007         }
1008 /*        if (ii[i] < nset)
1009             fprintf(fp,"Found %s data.\n",my_ener[i]);
1010  */ }
1011     /* Compute it all! */
1012     vvhh = alpha = kappa = cp = dcp = cv = NOTSET;
1013
1014     /* Temperature */
1015     tt = NOTSET;
1016     if (ii[eTemp] < nset)
1017     {
1018         tt    = edat->s[ii[eTemp]].av;
1019     }
1020     /* Volume */
1021     vv = varv = NOTSET;
1022     if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1023     {
1024         vv    = edat->s[ii[eVol]].av*NANO3;
1025         varv  = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1026         kappa = (varv/vv)/(BOLTZMANN*tt);
1027     }
1028     /* Enthalpy */
1029     hh = varh = NOTSET;
1030     if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1031     {
1032         hh    = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1033         varh  = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1034         cp    = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1035     }
1036     /* Total energy */
1037     et = varet = NOTSET;
1038     if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1039     {
1040         /* Only compute cv in constant volume runs, which we can test
1041            by checking whether the enthalpy was computed.
1042          */
1043         et    = edat->s[ii[eEtot]].av;
1044         varet = sqr(edat->s[ii[eEtot]].rmsd);
1045         cv    = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1046     }
1047     /* Alpha, dcp */
1048     if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1049     {
1050         vvhh = 0;
1051         for (j = 0; (j < edat->nframes); j++)
1052         {
1053             v     = edat->s[ii[eVol]].ener[j]*NANO3;
1054             h     = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
1055             vvhh += (v*h);
1056         }
1057         vvhh /= edat->nframes;
1058         alpha = (vvhh-vv*hh)/(vv*BOLTZMANN*tt*tt);
1059         dcp   = (vv*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1060     }
1061
1062     if (tt != NOTSET)
1063     {
1064         if (nmol < 2)
1065         {
1066             fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1067                     nmol);
1068         }
1069         fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1070         fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1071         fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1072         fprintf(fp, "please use the g_dos program.\n\n");
1073         fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1074                 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1075
1076         if (debug != NULL)
1077         {
1078             if (varv != NOTSET)
1079             {
1080                 fprintf(fp, "varv  =  %10g (m^6)\n", varv*AVOGADRO/nmol);
1081             }
1082             if (vvhh != NOTSET)
1083             {
1084                 fprintf(fp, "vvhh  =  %10g (m^3 J)\n", vvhh);
1085             }
1086         }
1087         if (vv != NOTSET)
1088         {
1089             fprintf(fp, "Volume                                   = %10g m^3/mol\n",
1090                     vv*AVOGADRO/nmol);
1091         }
1092         if (varh != NOTSET)
1093         {
1094             fprintf(fp, "Enthalpy                                 = %10g kJ/mol\n",
1095                     hh*AVOGADRO/(KILO*nmol));
1096         }
1097         if (alpha != NOTSET)
1098         {
1099             fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
1100                     alpha);
1101         }
1102         if (kappa != NOTSET)
1103         {
1104             fprintf(fp, "Isothermal Compressibility Kappa         = %10g (J/m^3)\n",
1105                     kappa);
1106             fprintf(fp, "Adiabatic bulk modulus                   = %10g (m^3/J)\n",
1107                     1.0/kappa);
1108         }
1109         if (cp != NOTSET)
1110         {
1111             fprintf(fp, "Heat capacity at constant pressure Cp    = %10g J/mol K\n",
1112                     cp);
1113         }
1114         if (cv != NOTSET)
1115         {
1116             fprintf(fp, "Heat capacity at constant volume Cv      = %10g J/mol K\n",
1117                     cv);
1118         }
1119         if (dcp != NOTSET)
1120         {
1121             fprintf(fp, "Cp-Cv                                    =  %10g J/mol K\n",
1122                     dcp);
1123         }
1124         please_cite(fp, "Allen1987a");
1125     }
1126     else
1127     {
1128         fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
1129     }
1130 }
1131
1132 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
1133                          gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
1134                          gmx_bool bVisco, const char *visfn, int nmol,
1135                          gmx_int64_t start_step, double start_t,
1136                          gmx_int64_t step, double t,
1137                          double time[], real reftemp,
1138                          enerdata_t *edat,
1139                          int nset, int set[], gmx_bool *bIsEner,
1140                          char **leg, gmx_enxnm_t *enm,
1141                          real Vaver, real ezero,
1142                          int nbmin, int nbmax,
1143                          const output_env_t oenv)
1144 {
1145     FILE           *fp;
1146     /* Check out the printed manual for equations! */
1147     double          Dt, aver, stddev, errest, delta_t, totaldrift;
1148     enerdata_t     *esum = NULL;
1149     real            xxx, integral, intBulk, Temp = 0, Pres = 0;
1150     real            sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
1151     double          beta = 0, expE, expEtot, *fee = NULL;
1152     gmx_int64_t     nsteps;
1153     int             nexact, nnotexact;
1154     double          x1m, x1mk;
1155     int             i, j, k, nout;
1156     real            chi2;
1157     char            buf[256], eebuf[100];
1158
1159     nsteps  = step - start_step + 1;
1160     if (nsteps < 1)
1161     {
1162         fprintf(stdout, "Not enough steps (%s) for statistics\n",
1163                 gmx_step_str(nsteps, buf));
1164     }
1165     else
1166     {
1167         /* Calculate the time difference */
1168         delta_t = t - start_t;
1169
1170         fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1171                 gmx_step_str(nsteps, buf), start_t, t, nset);
1172
1173         calc_averages(nset, edat, nbmin, nbmax);
1174
1175         if (bSum)
1176         {
1177             esum = calc_sum(nset, edat, nbmin, nbmax);
1178         }
1179
1180         if (edat->npoints == 0)
1181         {
1182             nexact    = 0;
1183             nnotexact = nset;
1184         }
1185         else
1186         {
1187             nexact    = 0;
1188             nnotexact = 0;
1189             for (i = 0; (i < nset); i++)
1190             {
1191                 if (edat->s[i].bExactStat)
1192                 {
1193                     nexact++;
1194                 }
1195                 else
1196                 {
1197                     nnotexact++;
1198                 }
1199             }
1200         }
1201
1202         if (nnotexact == 0)
1203         {
1204             fprintf(stdout, "All statistics are over %s points\n",
1205                     gmx_step_str(edat->npoints, buf));
1206         }
1207         else if (nexact == 0 || edat->npoints == edat->nframes)
1208         {
1209             fprintf(stdout, "All statistics are over %d points (frames)\n",
1210                     edat->nframes);
1211         }
1212         else
1213         {
1214             fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
1215             for (i = 0; (i < nset); i++)
1216             {
1217                 if (!edat->s[i].bExactStat)
1218                 {
1219                     fprintf(stdout, " '%s'", leg[i]);
1220                 }
1221             }
1222             fprintf(stdout, " %s has statistics over %d points (frames)\n",
1223                     nnotexact == 1 ? "is" : "are", edat->nframes);
1224             fprintf(stdout, "All other statistics are over %s points\n",
1225                     gmx_step_str(edat->npoints, buf));
1226         }
1227         fprintf(stdout, "\n");
1228
1229         fprintf(stdout, "%-24s %10s %10s %10s %10s",
1230                 "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
1231         if (bFee)
1232         {
1233             fprintf(stdout, "  %10s\n", "-kT ln<e^(E/kT)>");
1234         }
1235         else
1236         {
1237             fprintf(stdout, "\n");
1238         }
1239         fprintf(stdout, "-------------------------------------------------------------------------------\n");
1240
1241         /* Initiate locals, only used with -sum */
1242         expEtot = 0;
1243         if (bFee)
1244         {
1245             beta = 1.0/(BOLTZ*reftemp);
1246             snew(fee, nset);
1247         }
1248         for (i = 0; (i < nset); i++)
1249         {
1250             aver   = edat->s[i].av;
1251             stddev = edat->s[i].rmsd;
1252             errest = edat->s[i].ee;
1253
1254             if (bFee)
1255             {
1256                 expE = 0;
1257                 for (j = 0; (j < edat->nframes); j++)
1258                 {
1259                     expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1260                 }
1261                 if (bSum)
1262                 {
1263                     expEtot += expE/edat->nframes;
1264                 }
1265
1266                 fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
1267             }
1268             if (strstr(leg[i], "empera") != NULL)
1269             {
1270                 Temp = aver;
1271             }
1272             else if (strstr(leg[i], "olum") != NULL)
1273             {
1274                 Vaver = aver;
1275             }
1276             else if (strstr(leg[i], "essure") != NULL)
1277             {
1278                 Pres = aver;
1279             }
1280             if (bIsEner[i])
1281             {
1282                 pr_aver   = aver/nmol-ezero;
1283                 pr_stddev = stddev/nmol;
1284                 pr_errest = errest/nmol;
1285             }
1286             else
1287             {
1288                 pr_aver   = aver;
1289                 pr_stddev = stddev;
1290                 pr_errest = errest;
1291             }
1292
1293             /* Multiply the slope in steps with the number of steps taken */
1294             totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1295             if (bIsEner[i])
1296             {
1297                 totaldrift /= nmol;
1298             }
1299
1300             fprintf(stdout, "%-24s %10g %10s %10g %10g",
1301                     leg[i], pr_aver, ee_pr(pr_errest, eebuf), pr_stddev, totaldrift);
1302             if (bFee)
1303             {
1304                 fprintf(stdout, "  %10g", fee[i]);
1305             }
1306
1307             fprintf(stdout, "  (%s)\n", enm[set[i]].unit);
1308
1309             if (bFluct)
1310             {
1311                 for (j = 0; (j < edat->nframes); j++)
1312                 {
1313                     edat->s[i].ener[j] -= aver;
1314                 }
1315             }
1316         }
1317         if (bSum)
1318         {
1319             totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1320             fprintf(stdout, "%-24s %10g %10s %10s %10g  (%s)",
1321                     "Total", esum->s[0].av/nmol, ee_pr(esum->s[0].ee/nmol, eebuf),
1322                     "--", totaldrift/nmol, enm[set[0]].unit);
1323             /* pr_aver,pr_stddev,a,totaldrift */
1324             if (bFee)
1325             {
1326                 fprintf(stdout, "  %10g  %10g\n",
1327                         log(expEtot)/beta + esum->s[0].av/nmol, log(expEtot)/beta);
1328             }
1329             else
1330             {
1331                 fprintf(stdout, "\n");
1332             }
1333         }
1334
1335         /* Do correlation function */
1336         if (edat->nframes > 1)
1337         {
1338             Dt = delta_t/(edat->nframes - 1);
1339         }
1340         else
1341         {
1342             Dt = 0;
1343         }
1344         if (bVisco)
1345         {
1346             const char* leg[] = { "Shear", "Bulk" };
1347             real        factor;
1348             real      **eneset;
1349             real      **enesum;
1350
1351             /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1352
1353             /* Symmetrise tensor! (and store in first three elements)
1354              * And subtract average pressure!
1355              */
1356             snew(eneset, 12);
1357             for (i = 0; i < 12; i++)
1358             {
1359                 snew(eneset[i], edat->nframes);
1360             }
1361             snew(enesum, 3);
1362             for (i = 0; i < 3; i++)
1363             {
1364                 snew(enesum[i], edat->nframes);
1365             }
1366             for (i = 0; (i < edat->nframes); i++)
1367             {
1368                 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1369                 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1370                 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1371                 for (j = 3; j <= 11; j++)
1372                 {
1373                     eneset[j][i] = edat->s[j].ener[i];
1374                 }
1375                 eneset[11][i] -= Pres;
1376                 enesum[0][i]   = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
1377                 enesum[1][i]   = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
1378                 enesum[2][i]   = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
1379             }
1380
1381             einstein_visco("evisco.xvg", "eviscoi.xvg",
1382                            3, edat->nframes, enesum, Vaver, Temp, nsteps, time, oenv);
1383
1384             /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1385             /* Do it for shear viscosity */
1386             strcpy(buf, "Shear Viscosity");
1387             low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
1388                             (edat->nframes+1)/2, eneset, Dt,
1389                             eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1390
1391             /* Now for bulk viscosity */
1392             strcpy(buf, "Bulk Viscosity");
1393             low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1394                             (edat->nframes+1)/2, &(eneset[11]), Dt,
1395                             eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1396
1397             factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1398             fp     = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1399             xvgr_legend(fp, asize(leg), leg, oenv);
1400
1401             /* Use trapezium rule for integration */
1402             integral = 0;
1403             intBulk  = 0;
1404             nout     = get_acfnout();
1405             if ((nout < 2) || (nout >= edat->nframes/2))
1406             {
1407                 nout = edat->nframes/2;
1408             }
1409             for (i = 1; (i < nout); i++)
1410             {
1411                 integral += 0.5*(eneset[0][i-1]  + eneset[0][i])*factor;
1412                 intBulk  += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1413                 fprintf(fp, "%10g  %10g  %10g\n", (i*Dt), integral, intBulk);
1414             }
1415             ffclose(fp);
1416         }
1417         else if (bCorr)
1418         {
1419             if (bFluct)
1420             {
1421                 strcpy(buf, "Autocorrelation of Energy Fluctuations");
1422             }
1423             else
1424             {
1425                 strcpy(buf, "Energy Autocorrelation");
1426             }
1427 #if 0
1428             do_autocorr(corrfn, oenv, buf, edat->nframes,
1429                         bSum ? 1                 : nset,
1430                         bSum ? &edat->s[nset-1].ener : eneset,
1431                         (delta_t/edat->nframes), eacNormal, FALSE);
1432 #endif
1433         }
1434     }
1435 }
1436
1437 static void print_time(FILE *fp, double t)
1438 {
1439     fprintf(fp, "%12.6f", t);
1440 }
1441
1442 static void print1(FILE *fp, gmx_bool bDp, real e)
1443 {
1444     if (bDp)
1445     {
1446         fprintf(fp, "  %16.12f", e);
1447     }
1448     else
1449     {
1450         fprintf(fp, "  %10.6f", e);
1451     }
1452 }
1453
1454 static void fec(const char *ene2fn, const char *runavgfn,
1455                 real reftemp, int nset, int set[], char *leg[],
1456                 enerdata_t *edat, double time[],
1457                 const output_env_t oenv)
1458 {
1459     const char * ravgleg[] = {
1460         "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1461         "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
1462     };
1463     FILE        *fp;
1464     ener_file_t  enx;
1465     int          nre, timecheck, step, nenergy, nenergy2, maxenergy;
1466     int          i, j;
1467     gmx_bool     bCont;
1468     real         aver, beta;
1469     real       **eneset2;
1470     double       dE, sum;
1471     gmx_enxnm_t *enm = NULL;
1472     t_enxframe  *fr;
1473     char         buf[22];
1474
1475     /* read second energy file */
1476     snew(fr, 1);
1477     enm = NULL;
1478     enx = open_enx(ene2fn, "r");
1479     do_enxnms(enx, &(fr->nre), &enm);
1480
1481     snew(eneset2, nset+1);
1482     nenergy2  = 0;
1483     maxenergy = 0;
1484     timecheck = 0;
1485     do
1486     {
1487         /* This loop searches for the first frame (when -b option is given),
1488          * or when this has been found it reads just one energy frame
1489          */
1490         do
1491         {
1492             bCont = do_enx(enx, fr);
1493
1494             if (bCont)
1495             {
1496                 timecheck = check_times(fr->t);
1497             }
1498
1499         }
1500         while (bCont && (timecheck < 0));
1501
1502         /* Store energies for analysis afterwards... */
1503         if ((timecheck == 0) && bCont)
1504         {
1505             if (fr->nre > 0)
1506             {
1507                 if (nenergy2 >= maxenergy)
1508                 {
1509                     maxenergy += 1000;
1510                     for (i = 0; i <= nset; i++)
1511                     {
1512                         srenew(eneset2[i], maxenergy);
1513                     }
1514                 }
1515                 if (fr->t != time[nenergy2])
1516                 {
1517                     fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
1518                             fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
1519                 }
1520                 for (i = 0; i < nset; i++)
1521                 {
1522                     eneset2[i][nenergy2] = fr->ener[set[i]].e;
1523                 }
1524                 nenergy2++;
1525             }
1526         }
1527     }
1528     while (bCont && (timecheck == 0));
1529
1530     /* check */
1531     if (edat->nframes != nenergy2)
1532     {
1533         fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
1534                 edat->nframes, nenergy2);
1535     }
1536     nenergy = min(edat->nframes, nenergy2);
1537
1538     /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1539     fp = NULL;
1540     if (runavgfn)
1541     {
1542         fp = xvgropen(runavgfn, "Running average free energy difference",
1543                       "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
1544         xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1545     }
1546     fprintf(stdout, "\n%-24s %10s\n",
1547             "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1548     sum  = 0;
1549     beta = 1.0/(BOLTZ*reftemp);
1550     for (i = 0; i < nset; i++)
1551     {
1552         if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1553         {
1554             fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
1555                     leg[i], enm[set[i]].name);
1556         }
1557         for (j = 0; j < nenergy; j++)
1558         {
1559             dE   = eneset2[i][j] - edat->s[i].ener[j];
1560             sum += exp(-dE*beta);
1561             if (fp)
1562             {
1563                 fprintf(fp, "%10g %10g %10g\n",
1564                         time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1565             }
1566         }
1567         aver = -BOLTZ*reftemp*log(sum/nenergy);
1568         fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1569     }
1570     if (fp)
1571     {
1572         ffclose(fp);
1573     }
1574     sfree(fr);
1575 }
1576
1577
1578 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl,
1579                     const char *filename, gmx_bool bDp,
1580                     int *blocks, int *hists, int *samples, int *nlambdas,
1581                     const output_env_t oenv)
1582 {
1583     const char  *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1584     char         title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1585     char         buf[STRLEN];
1586     gmx_bool     first       = FALSE;
1587     int          nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1588     int          i, j, k;
1589     /* coll data */
1590     double       temp              = 0, start_time = 0, delta_time = 0, start_lambda = 0, delta_lambda = 0;
1591     static int   setnr             = 0;
1592     double      *native_lambda_vec = NULL;
1593     const char **lambda_components = NULL;
1594     int          n_lambda_vec      = 0;
1595     gmx_bool     changing_lambda   = FALSE;
1596     int          lambda_fep_state;
1597
1598     /* now count the blocks & handle the global dh data */
1599     for (i = 0; i < fr->nblock; i++)
1600     {
1601         if (fr->block[i].id == enxDHHIST)
1602         {
1603             nblock_hist++;
1604         }
1605         else if (fr->block[i].id == enxDH)
1606         {
1607             nblock_dh++;
1608         }
1609         else if (fr->block[i].id == enxDHCOLL)
1610         {
1611             nblock_dhcoll++;
1612             if ( (fr->block[i].nsub < 1) ||
1613                  (fr->block[i].sub[0].type != xdr_datatype_double) ||
1614                  (fr->block[i].sub[0].nr < 5))
1615             {
1616                 gmx_fatal(FARGS, "Unexpected block data");
1617             }
1618
1619             /* read the data from the DHCOLL block */
1620             temp            =         fr->block[i].sub[0].dval[0];
1621             start_time      =   fr->block[i].sub[0].dval[1];
1622             delta_time      =   fr->block[i].sub[0].dval[2];
1623             start_lambda    = fr->block[i].sub[0].dval[3];
1624             delta_lambda    = fr->block[i].sub[0].dval[4];
1625             changing_lambda = (delta_lambda != 0);
1626             if (fr->block[i].nsub > 1)
1627             {
1628                 lambda_fep_state = fr->block[i].sub[1].ival[0];
1629                 if (n_lambda_vec == 0)
1630                 {
1631                     n_lambda_vec = fr->block[i].sub[1].ival[1];
1632                 }
1633                 else
1634                 {
1635                     if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1636                     {
1637                         gmx_fatal(FARGS,
1638                                   "Unexpected change of basis set in lambda");
1639                     }
1640                 }
1641                 if (lambda_components == NULL)
1642                 {
1643                     snew(lambda_components, n_lambda_vec);
1644                 }
1645                 if (native_lambda_vec == NULL)
1646                 {
1647                     snew(native_lambda_vec, n_lambda_vec);
1648                 }
1649                 for (j = 0; j < n_lambda_vec; j++)
1650                 {
1651                     native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1652                     lambda_components[j] =
1653                         efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1654                 }
1655             }
1656         }
1657     }
1658
1659     if (nblock_hist == 0 && nblock_dh == 0)
1660     {
1661         /* don't do anything */
1662         return;
1663     }
1664     if (nblock_hist > 0 && nblock_dh > 0)
1665     {
1666         gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1667     }
1668     if (!*fp_dhdl)
1669     {
1670         if (nblock_dh > 0)
1671         {
1672             /* we have standard, non-histogram data --
1673                call open_dhdl to open the file */
1674             /* TODO this is an ugly hack that needs to be fixed: this will only
1675                work if the order of data is always the same and if we're
1676                only using the g_energy compiled with the mdrun that produced
1677                the ener.edr. */
1678             *fp_dhdl = open_dhdl(filename, ir, oenv);
1679         }
1680         else
1681         {
1682             sprintf(title, "N(%s)", deltag);
1683             sprintf(label_x, "%s (%s)", deltag, unit_energy);
1684             sprintf(label_y, "Samples");
1685             *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1686             sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1687             xvgr_subtitle(*fp_dhdl, buf, oenv);
1688         }
1689     }
1690
1691     (*hists)   += nblock_hist;
1692     (*blocks)  += nblock_dh;
1693     (*nlambdas) = nblock_hist+nblock_dh;
1694
1695     /* write the data */
1696     if (nblock_hist > 0)
1697     {
1698         gmx_int64_t sum = 0;
1699         /* histograms */
1700         for (i = 0; i < fr->nblock; i++)
1701         {
1702             t_enxblock *blk = &(fr->block[i]);
1703             if (blk->id == enxDHHIST)
1704             {
1705                 double          foreign_lambda, dx;
1706                 gmx_int64_t     x0;
1707                 int             nhist, derivative;
1708
1709                 /* check the block types etc. */
1710                 if ( (blk->nsub < 2) ||
1711                      (blk->sub[0].type != xdr_datatype_double) ||
1712                      (blk->sub[1].type != xdr_datatype_large_int) ||
1713                      (blk->sub[0].nr < 2)  ||
1714                      (blk->sub[1].nr < 2) )
1715                 {
1716                     gmx_fatal(FARGS, "Unexpected block data in file");
1717                 }
1718                 foreign_lambda = blk->sub[0].dval[0];
1719                 dx             = blk->sub[0].dval[1];
1720                 nhist          = blk->sub[1].lval[0];
1721                 derivative     = blk->sub[1].lval[1];
1722                 for (j = 0; j < nhist; j++)
1723                 {
1724                     const char *lg[1];
1725                     x0 = blk->sub[1].lval[2+j];
1726
1727                     if (!derivative)
1728                     {
1729                         sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1730                                 deltag, lambda, foreign_lambda,
1731                                 lambda, start_lambda);
1732                     }
1733                     else
1734                     {
1735                         sprintf(legend, "N(%s | %s=%g)",
1736                                 dhdl, lambda, start_lambda);
1737                     }
1738
1739                     lg[0] = legend;
1740                     xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1741                     setnr++;
1742                     for (k = 0; k < blk->sub[j+2].nr; k++)
1743                     {
1744                         int    hist;
1745                         double xmin, xmax;
1746
1747                         hist = blk->sub[j+2].ival[k];
1748                         xmin = (x0+k)*dx;
1749                         xmax = (x0+k+1)*dx;
1750                         fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1751                                 xmax, hist);
1752                         sum += hist;
1753                     }
1754                     /* multiple histogram data blocks in one histogram
1755                        mean that the second one is the reverse of the first one:
1756                        for dhdl derivatives, it's important to know both the
1757                        maximum and minimum values */
1758                     dx = -dx;
1759                 }
1760             }
1761         }
1762         (*samples) += (int)(sum/nblock_hist);
1763     }
1764     else
1765     {
1766         /* raw dh */
1767         int    len      = 0;
1768         char **setnames = NULL;
1769         int    nnames   = nblock_dh;
1770
1771         for (i = 0; i < fr->nblock; i++)
1772         {
1773             t_enxblock *blk = &(fr->block[i]);
1774             if (blk->id == enxDH)
1775             {
1776                 if (len == 0)
1777                 {
1778                     len = blk->sub[2].nr;
1779                 }
1780                 else
1781                 {
1782                     if (len != blk->sub[2].nr)
1783                     {
1784                         gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1785                     }
1786                 }
1787             }
1788         }
1789         (*samples) += len;
1790
1791         for (i = 0; i < len; i++)
1792         {
1793             double time = start_time + delta_time*i;
1794
1795             fprintf(*fp_dhdl, "%.4f ", time);
1796
1797             for (j = 0; j < fr->nblock; j++)
1798             {
1799                 t_enxblock *blk = &(fr->block[j]);
1800                 if (blk->id == enxDH)
1801                 {
1802                     double value;
1803                     if (blk->sub[2].type == xdr_datatype_float)
1804                     {
1805                         value = blk->sub[2].fval[i];
1806                     }
1807                     else
1808                     {
1809                         value = blk->sub[2].dval[i];
1810                     }
1811                     /* we need to decide which data type it is based on the count*/
1812
1813                     if (j == 1 && ir->bExpanded)
1814                     {
1815                         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! */
1816                     }
1817                     else
1818                     {
1819                         if (bDp)
1820                         {
1821                             fprintf(*fp_dhdl, " %#.12g", value);   /* print normal precision */
1822                         }
1823                         else
1824                         {
1825                             fprintf(*fp_dhdl, " %#.8g", value);   /* print normal precision */
1826                         }
1827                     }
1828                 }
1829             }
1830             fprintf(*fp_dhdl, "\n");
1831         }
1832     }
1833 }
1834
1835
1836 int gmx_energy(int argc, char *argv[])
1837 {
1838     const char        *desc[] = {
1839         "[THISMODULE] 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 [THISMODULE] 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 [gmx-dos] 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_int64_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 }