src/tools/gmx_energy.c

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