src/gromacs/gmxana/gmx_energy.c

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