src/gromacs/gmxana/gmx_energy.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include <math.h>
  42 #include <string.h>
  43
  44 #include "typedefs.h"
  45 #include "gmx_fatal.h"
  46 #include "vec.h"
  47 #include "string2.h"
  48 #include "gromacs/utility/smalloc.h"
  49 #include "gromacs/fileio/enxio.h"
  50 #include "gromacs/commandline/pargs.h"
  51 #include "names.h"
  52 #include "copyrite.h"
  53 #include "macros.h"
  54 #include "xvgr.h"
  55 #include "gstat.h"
  56 #include "physics.h"
  57 #include "gromacs/fileio/tpxio.h"
  58 #include "gromacs/fileio/trxio.h"
  59 #include "viewit.h"
  60 #include "mtop_util.h"
  61 #include "gmx_ana.h"
  62 #include "mdebin.h"
  63
  64 static real       minthird = -1.0/3.0, minsixth = -1.0/6.0;
  65
  66 typedef struct {
  67     real sum;
  68     real sum2;
  69 } exactsum_t;
  70
  71 typedef struct {
  72     real       *ener;
  73     exactsum_t *es;
  74     gmx_bool    bExactStat;
  75     double      av;
  76     double      rmsd;
  77     double      ee;
  78     double      slope;
  79 } enerdat_t;
  80
  81 typedef struct {
  82     gmx_int64_t      nsteps;
  83     gmx_int64_t      npoints;
  84     int              nframes;
  85     int             *step;
  86     int             *steps;
  87     int             *points;
  88     enerdat_t       *s;
  89 } enerdata_t;
  90
  91 static double mypow(double x, double y)
  92 {
  93     if (x > 0)
  94     {
  95         return pow(x, y);
  96     }
  97     else
  98     {
  99         return 0.0;
 100     }
 101 }
 102
 103 static int *select_it(int nre, char *nm[], int *nset)
 104 {
 105     gmx_bool *bE;
 106     int       n, k, j, i;
 107     int      *set;
 108     gmx_bool  bVerbose = TRUE;
 109
 110     if ((getenv("VERBOSE")) != NULL)
 111     {
 112         bVerbose = FALSE;
 113     }
 114
 115     fprintf(stderr, "Select the terms you want from the following list\n");
 116     fprintf(stderr, "End your selection with 0\n");
 117
 118     if (bVerbose)
 119     {
 120         for (k = 0; (k < nre); )
 121         {
 122             for (j = 0; (j < 4) && (k < nre); j++, k++)
 123             {
 124                 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
 125             }
 126             fprintf(stderr, "\n");
 127         }
 128     }
 129
 130     snew(bE, nre);
 131     do
 132     {
 133         if (1 != scanf("%d", &n))
 134         {
 135             gmx_fatal(FARGS, "Error reading user input");
 136         }
 137         if ((n > 0) && (n <= nre))
 138         {
 139             bE[n-1] = TRUE;
 140         }
 141     }
 142     while (n != 0);
 143
 144     snew(set, nre);
 145     for (i = (*nset) = 0; (i < nre); i++)
 146     {
 147         if (bE[i])
 148         {
 149             set[(*nset)++] = i;
 150         }
 151     }
 152
 153     sfree(bE);
 154
 155     return set;
 156 }
 157
 158 static void chomp(char *buf)
 159 {
 160     int len = strlen(buf);
 161
 162     while ((len > 0) && (buf[len-1] == '\n'))
 163     {
 164         buf[len-1] = '\0';
 165         len--;
 166     }
 167 }
 168
 169 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
 170 {
 171     gmx_bool   *bE;
 172     int         n, k, kk, j, i, nmatch, nind, nss;
 173     int        *set;
 174     gmx_bool    bEOF, bVerbose = TRUE, bLong = FALSE;
 175     char       *ptr, buf[STRLEN];
 176     const char *fm4   = "%3d  %-14s";
 177     const char *fm2   = "%3d  %-34s";
 178     char      **newnm = NULL;
 179
 180     if ((getenv("VERBOSE")) != NULL)
 181     {
 182         bVerbose = FALSE;
 183     }
 184
 185     fprintf(stderr, "\n");
 186     fprintf(stderr, "Select the terms you want from the following list by\n");
 187     fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
 188     fprintf(stderr, "End your selection with an empty line or a zero.\n");
 189     fprintf(stderr, "-------------------------------------------------------------------\n");
 190
 191     snew(newnm, nre);
 192     j = 0;
 193     for (k = 0; k < nre; k++)
 194     {
 195         newnm[k] = strdup(nm[k].name);
 196         /* Insert dashes in all the names */
 197         while ((ptr = strchr(newnm[k], ' ')) != NULL)
 198         {
 199             *ptr = '-';
 200         }
 201         if (bVerbose)
 202         {
 203             if (j == 0)
 204             {
 205                 if (k > 0)
 206                 {
 207                     fprintf(stderr, "\n");
 208                 }
 209                 bLong = FALSE;
 210                 for (kk = k; kk < k+4; kk++)
 211                 {
 212                     if (kk < nre && strlen(nm[kk].name) > 14)
 213                     {
 214                         bLong = TRUE;
 215                     }
 216                 }
 217             }
 218             else
 219             {
 220                 fprintf(stderr, " ");
 221             }
 222             if (!bLong)
 223             {
 224                 fprintf(stderr, fm4, k+1, newnm[k]);
 225                 j++;
 226                 if (j == 4)
 227                 {
 228                     j = 0;
 229                 }
 230             }
 231             else
 232             {
 233                 fprintf(stderr, fm2, k+1, newnm[k]);
 234                 j++;
 235                 if (j == 2)
 236                 {
 237                     j = 0;
 238                 }
 239             }
 240         }
 241     }
 242     if (bVerbose)
 243     {
 244         fprintf(stderr, "\n\n");
 245     }
 246
 247     snew(bE, nre);
 248
 249     bEOF = FALSE;
 250     while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
 251     {
 252         /* Remove newlines */
 253         chomp(buf);
 254
 255         /* Remove spaces */
 256         trim(buf);
 257
 258         /* Empty line means end of input */
 259         bEOF = (strlen(buf) == 0);
 260         if (!bEOF)
 261         {
 262             ptr = buf;
 263             do
 264             {
 265                 if (!bEOF)
 266                 {
 267                     /* First try to read an integer */
 268                     nss   = sscanf(ptr, "%d", &nind);
 269                     if (nss == 1)
 270                     {
 271                         /* Zero means end of input */
 272                         if (nind == 0)
 273                         {
 274                             bEOF = TRUE;
 275                         }
 276                         else if ((1 <= nind) && (nind <= nre))
 277                         {
 278                             bE[nind-1] = TRUE;
 279                         }
 280                         else
 281                         {
 282                             fprintf(stderr, "number %d is out of range\n", nind);
 283                         }
 284                     }
 285                     else
 286                     {
 287                         /* Now try to read a string */
 288                         i      = strlen(ptr);
 289                         nmatch = 0;
 290                         for (nind = 0; nind < nre; nind++)
 291                         {
 292                             if (gmx_strcasecmp(newnm[nind], ptr) == 0)
 293                             {
 294                                 bE[nind] = TRUE;
 295                                 nmatch++;
 296                             }
 297                         }
 298                         if (nmatch == 0)
 299                         {
 300                             i      = strlen(ptr);
 301                             nmatch = 0;
 302                             for (nind = 0; nind < nre; nind++)
 303                             {
 304                                 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
 305                                 {
 306                                     bE[nind] = TRUE;
 307                                     nmatch++;
 308                                 }
 309                             }
 310                             if (nmatch == 0)
 311                             {
 312                                 fprintf(stderr, "String '%s' does not match anything\n", ptr);
 313                             }
 314                         }
 315                     }
 316                 }
 317                 /* Look for the first space, and remove spaces from there */
 318                 if ((ptr = strchr(ptr, ' ')) != NULL)
 319                 {
 320                     trim(ptr);
 321                 }
 322             }
 323             while (!bEOF && (ptr && (strlen(ptr) > 0)));
 324         }
 325     }
 326
 327     snew(set, nre);
 328     for (i = (*nset) = 0; (i < nre); i++)
 329     {
 330         if (bE[i])
 331         {
 332             set[(*nset)++] = i;
 333         }
 334     }
 335
 336     sfree(bE);
 337
 338     if (*nset == 0)
 339     {
 340         gmx_fatal(FARGS, "No energy terms selected");
 341     }
 342
 343     for (i = 0; (i < nre); i++)
 344     {
 345         sfree(newnm[i]);
 346     }
 347     sfree(newnm);
 348
 349     return set;
 350 }
 351
 352 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
 353 {
 354     gmx_mtop_t  mtop;
 355     int         natoms;
 356     t_iatom    *iatom;
 357     matrix      box;
 358
 359     /* all we need is the ir to be able to write the label */
 360     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
 361 }
 362
 363 static void get_orires_parms(const char *topnm,
 364                              int *nor, int *nex, int **label, real **obs)
 365 {
 366     gmx_mtop_t      mtop;
 367     gmx_localtop_t *top;
 368     t_inputrec      ir;
 369     t_iparams      *ip;
 370     int             natoms, i;
 371     t_iatom        *iatom;
 372     int             nb;
 373     matrix          box;
 374
 375     read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
 376     top = gmx_mtop_generate_local_top(&mtop, &ir);
 377
 378     ip       = top->idef.iparams;
 379     iatom    = top->idef.il[F_ORIRES].iatoms;
 380
 381     /* Count how many distance restraint there are... */
 382     nb = top->idef.il[F_ORIRES].nr;
 383     if (nb == 0)
 384     {
 385         gmx_fatal(FARGS, "No orientation restraints in topology!\n");
 386     }
 387
 388     *nor = nb/3;
 389     *nex = 0;
 390     snew(*label, *nor);
 391     snew(*obs, *nor);
 392     for (i = 0; i < nb; i += 3)
 393     {
 394         (*label)[i/3] = ip[iatom[i]].orires.label;
 395         (*obs)[i/3]   = ip[iatom[i]].orires.obs;
 396         if (ip[iatom[i]].orires.ex >= *nex)
 397         {
 398             *nex = ip[iatom[i]].orires.ex+1;
 399         }
 400     }
 401     fprintf(stderr, "Found %d orientation restraints with %d experiments",
 402             *nor, *nex);
 403 }
 404
 405 static int get_bounds(const char *topnm,
 406                       real **bounds, int **index, int **dr_pair, int *npairs,
 407                       gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
 408 {
 409     gmx_localtop_t *top;
 410     t_functype     *functype;
 411     t_iparams      *ip;
 412     int             natoms, i, j, k, type, ftype, natom;
 413     t_ilist        *disres;
 414     t_iatom        *iatom;
 415     real           *b;
 416     int            *ind, *pair;
 417     int             nb, label1;
 418     matrix          box;
 419
 420     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
 421     snew(*ltop, 1);
 422     top   = gmx_mtop_generate_local_top(mtop, ir);
 423     *ltop = top;
 424
 425     functype = top->idef.functype;
 426     ip       = top->idef.iparams;
 427
 428     /* Count how many distance restraint there are... */
 429     nb = top->idef.il[F_DISRES].nr;
 430     if (nb == 0)
 431     {
 432         gmx_fatal(FARGS, "No distance restraints in topology!\n");
 433     }
 434
 435     /* Allocate memory */
 436     snew(b, nb);
 437     snew(ind, nb);
 438     snew(pair, nb+1);
 439
 440     /* Fill the bound array */
 441     nb = 0;
 442     for (i = 0; (i < top->idef.ntypes); i++)
 443     {
 444         ftype = functype[i];
 445         if (ftype == F_DISRES)
 446         {
 447
 448             label1  = ip[i].disres.label;
 449             b[nb]   = ip[i].disres.up1;
 450             ind[nb] = label1;
 451             nb++;
 452         }
 453     }
 454     *bounds = b;
 455
 456     /* Fill the index array */
 457     label1  = -1;
 458     disres  = &(top->idef.il[F_DISRES]);
 459     iatom   = disres->iatoms;
 460     for (i = j = k = 0; (i < disres->nr); )
 461     {
 462         type  = iatom[i];
 463         ftype = top->idef.functype[type];
 464         natom = interaction_function[ftype].nratoms+1;
 465         if (label1 != top->idef.iparams[type].disres.label)
 466         {
 467             pair[j] = k;
 468             label1  = top->idef.iparams[type].disres.label;
 469             j++;
 470         }
 471         k++;
 472         i += natom;
 473     }
 474     pair[j]  = k;
 475     *npairs  = k;
 476     if (j != nb)
 477     {
 478         gmx_incons("get_bounds for distance restraints");
 479     }
 480
 481     *index   = ind;
 482     *dr_pair = pair;
 483
 484     return nb;
 485 }
 486
 487 static void calc_violations(real rt[], real rav3[], int nb, int index[],
 488                             real bounds[], real *viol, double *st, double *sa)
 489 {
 490     const   real sixth = 1.0/6.0;
 491     int          i, j;
 492     double       rsum, rav, sumaver, sumt;
 493
 494     sumaver = 0;
 495     sumt    = 0;
 496     for (i = 0; (i < nb); i++)
 497     {
 498         rsum = 0.0;
 499         rav  = 0.0;
 500         for (j = index[i]; (j < index[i+1]); j++)
 501         {
 502             if (viol)
 503             {
 504                 viol[j] += mypow(rt[j], -3.0);
 505             }
 506             rav     += sqr(rav3[j]);
 507             rsum    += mypow(rt[j], -6);
 508         }
 509         rsum    = max(0.0, mypow(rsum, -sixth)-bounds[i]);
 510         rav     = max(0.0, mypow(rav, -sixth)-bounds[i]);
 511
 512         sumt    += rsum;
 513         sumaver += rav;
 514     }
 515     *st = sumt;
 516     *sa = sumaver;
 517 }
 518
 519 static void analyse_disre(const char *voutfn,    int nframes,
 520                           real violaver[], real bounds[], int index[],
 521                           int pair[],      int nbounds,
 522                           const output_env_t oenv)
 523 {
 524     FILE   *vout;
 525     double  sum, sumt, sumaver;
 526     int     i, j;
 527
 528     /* Subtract bounds from distances, to calculate violations */
 529     calc_violations(violaver, violaver,
 530                     nbounds, pair, bounds, NULL, &sumt, &sumaver);
 531
 532 #ifdef DEBUG
 533     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 534             sumaver);
 535     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
 536             sumt);
 537 #endif
 538     vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
 539                     oenv);
 540     sum  = 0.0;
 541     sumt = 0.0;
 542     for (i = 0; (i < nbounds); i++)
 543     {
 544         /* Do ensemble averaging */
 545         sumaver = 0;
 546         for (j = pair[i]; (j < pair[i+1]); j++)
 547         {
 548             sumaver += sqr(violaver[j]/nframes);
 549         }
 550         sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
 551
 552         sumt   += sumaver;
 553         sum     = max(sum, sumaver);
 554         fprintf(vout, "%10d  %10.5e\n", index[i], sumaver);
 555     }
 556 #ifdef DEBUG
 557     for (j = 0; (j < dr.ndr); j++)
 558     {
 559         fprintf(vout, "%10d  %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
 560     }
 561 #endif
 562     gmx_ffclose(vout);
 563
 564     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 565             sumt);
 566     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
 567
 568     do_view(oenv, voutfn, "-graphtype bar");
 569 }
 570
 571 static void einstein_visco(const char *fn, const char *fni, int nsets,
 572                            int nframes, real **sum,
 573                            real V, real T, int nsteps, double time[],
 574                            const output_env_t oenv)
 575 {
 576     FILE  *fp0, *fp1;
 577     double av[4], avold[4];
 578     double fac, dt, di;
 579     int    i, j, m, nf4;
 580
 581     if (nframes < 1)
 582     {
 583         return;
 584     }
 585
 586     dt  = (time[1]-time[0]);
 587     nf4 = nframes/4+1;
 588
 589     for (i = 0; i <= nsets; i++)
 590     {
 591         avold[i] = 0;
 592     }
 593     fp0 = xvgropen(fni, "Shear viscosity integral",
 594                    "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
 595     fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
 596                    "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
 597     for (i = 1; i < nf4; i++)
 598     {
 599         fac = dt*nframes/nsteps;
 600         for (m = 0; m <= nsets; m++)
 601         {
 602             av[m] = 0;
 603         }
 604         for (j = 0; j < nframes-i; j++)
 605         {
 606             for (m = 0; m < nsets; m++)
 607             {
 608                 di   = sqr(fac*(sum[m][j+i]-sum[m][j]));
 609
 610                 av[m]     += di;
 611                 av[nsets] += di/nsets;
 612             }
 613         }
 614         /* Convert to SI for the viscosity */
 615         fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
 616         fprintf(fp0, "%10g", time[i]-time[0]);
 617         for (m = 0; (m <= nsets); m++)
 618         {
 619             av[m] = fac*av[m];
 620             fprintf(fp0, "  %10g", av[m]);
 621         }
 622         fprintf(fp0, "\n");
 623         fprintf(fp1, "%10g", 0.5*(time[i]+time[i-1])-time[0]);
 624         for (m = 0; (m <= nsets); m++)
 625         {
 626             fprintf(fp1, "  %10g", (av[m]-avold[m])/dt);
 627             avold[m] = av[m];
 628         }
 629         fprintf(fp1, "\n");
 630     }
 631     gmx_ffclose(fp0);
 632     gmx_ffclose(fp1);
 633 }
 634
 635 typedef struct {
 636     gmx_int64_t     np;
 637     double          sum;
 638     double          sav;
 639     double          sav2;
 640 } ee_sum_t;
 641
 642 typedef struct {
 643     int             b;
 644     ee_sum_t        sum;
 645     gmx_int64_t     nst;
 646     gmx_int64_t     nst_min;
 647 } ener_ee_t;
 648
 649 static void clear_ee_sum(ee_sum_t *ees)
 650 {
 651     ees->sav  = 0;
 652     ees->sav2 = 0;
 653     ees->np   = 0;
 654     ees->sum  = 0;
 655 }
 656
 657 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
 658 {
 659     ees->np  += np;
 660     ees->sum += sum;
 661 }
 662
 663 static void add_ee_av(ee_sum_t *ees)
 664 {
 665     double av;
 666
 667     av         = ees->sum/ees->np;
 668     ees->sav  += av;
 669     ees->sav2 += av*av;
 670     ees->np    = 0;
 671     ees->sum   = 0;
 672 }
 673
 674 static double calc_ee2(int nb, ee_sum_t *ees)
 675 {
 676     return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
 677 }
 678
 679 static void set_ee_av(ener_ee_t *eee)
 680 {
 681     if (debug)
 682     {
 683         char buf[STEPSTRSIZE];
 684         fprintf(debug, "Storing average for err.est.: %s steps\n",
 685                 gmx_step_str(eee->nst, buf));
 686     }
 687     add_ee_av(&eee->sum);
 688     eee->b++;
 689     if (eee->b == 1 || eee->nst < eee->nst_min)
 690     {
 691         eee->nst_min = eee->nst;
 692     }
 693     eee->nst = 0;
 694 }
 695
 696 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
 697 {
 698     int             nb, i, f, nee;
 699     double          sum, sum2, sump, see2;
 700     gmx_int64_t     steps, np, p, bound_nb;
 701     enerdat_t      *ed;
 702     exactsum_t     *es;
 703     gmx_bool        bAllZero;
 704     double          x, sx, sy, sxx, sxy;
 705     ener_ee_t      *eee;
 706
 707     /* Check if we have exact statistics over all points */
 708     for (i = 0; i < nset; i++)
 709     {
 710         ed             = &edat->s[i];
 711         ed->bExactStat = FALSE;
 712         if (edat->npoints > 0)
 713         {
 714             /* All energy file sum entries 0 signals no exact sums.
 715              * But if all energy values are 0, we still have exact sums.
 716              */
 717             bAllZero = TRUE;
 718             for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
 719             {
 720                 if (ed->ener[i] != 0)
 721                 {
 722                     bAllZero = FALSE;
 723                 }
 724                 ed->bExactStat = (ed->es[f].sum != 0);
 725             }
 726             if (bAllZero)
 727             {
 728                 ed->bExactStat = TRUE;
 729             }
 730         }
 731     }
 732
 733     snew(eee, nbmax+1);
 734     for (i = 0; i < nset; i++)
 735     {
 736         ed = &edat->s[i];
 737
 738         sum  = 0;
 739         sum2 = 0;
 740         np   = 0;
 741         sx   = 0;
 742         sy   = 0;
 743         sxx  = 0;
 744         sxy  = 0;
 745         for (nb = nbmin; nb <= nbmax; nb++)
 746         {
 747             eee[nb].b     = 0;
 748             clear_ee_sum(&eee[nb].sum);
 749             eee[nb].nst     = 0;
 750             eee[nb].nst_min = 0;
 751         }
 752         for (f = 0; f < edat->nframes; f++)
 753         {
 754             es = &ed->es[f];
 755
 756             if (ed->bExactStat)
 757             {
 758                 /* Add the sum and the sum of variances to the totals. */
 759                 p     = edat->points[f];
 760                 sump  = es->sum;
 761                 sum2 += es->sum2;
 762                 if (np > 0)
 763                 {
 764                     sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
 765                         *np*(np + p)/p;
 766                 }
 767             }
 768             else
 769             {
 770                 /* Add a single value to the sum and sum of squares. */
 771                 p     = 1;
 772                 sump  = ed->ener[f];
 773                 sum2 += dsqr(sump);
 774             }
 775
 776             /* sum has to be increased after sum2 */
 777             np  += p;
 778             sum += sump;
 779
 780             /* For the linear regression use variance 1/p.
 781              * Note that sump is the sum, not the average, so we don't need p*.
 782              */
 783             x    = edat->step[f] - 0.5*(edat->steps[f] - 1);
 784             sx  += p*x;
 785             sy  += sump;
 786             sxx += p*x*x;
 787             sxy += x*sump;
 788
 789             for (nb = nbmin; nb <= nbmax; nb++)
 790             {
 791                 /* Check if the current end step is closer to the desired
 792                  * block boundary than the next end step.
 793                  */
 794                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 795                 if (eee[nb].nst > 0 &&
 796                     bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
 797                 {
 798                     set_ee_av(&eee[nb]);
 799                 }
 800                 if (f == 0)
 801                 {
 802                     eee[nb].nst = 1;
 803                 }
 804                 else
 805                 {
 806                     eee[nb].nst += edat->step[f] - edat->step[f-1];
 807                 }
 808                 if (ed->bExactStat)
 809                 {
 810                     add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
 811                 }
 812                 else
 813                 {
 814                     add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
 815                 }
 816                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 817                 if (edat->step[f]*nb >= bound_nb)
 818                 {
 819                     set_ee_av(&eee[nb]);
 820                 }
 821             }
 822         }
 823
 824         edat->s[i].av = sum/np;
 825         if (ed->bExactStat)
 826         {
 827             edat->s[i].rmsd = sqrt(sum2/np);
 828         }
 829         else
 830         {
 831             edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
 832         }
 833
 834         if (edat->nframes > 1)
 835         {
 836             edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
 837         }
 838         else
 839         {
 840             edat->s[i].slope = 0;
 841         }
 842
 843         nee  = 0;
 844         see2 = 0;
 845         for (nb = nbmin; nb <= nbmax; nb++)
 846         {
 847             /* Check if we actually got nb blocks and if the smallest
 848              * block is not shorter than 80% of the average.
 849              */
 850             if (debug)
 851             {
 852                 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
 853                 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
 854                         nb, eee[nb].b,
 855                         gmx_step_str(eee[nb].nst_min, buf1),
 856                         gmx_step_str(edat->nsteps, buf2));
 857             }
 858             if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
 859             {
 860                 see2 += calc_ee2(nb, &eee[nb].sum);
 861                 nee++;
 862             }
 863         }
 864         if (nee > 0)
 865         {
 866             edat->s[i].ee = sqrt(see2/nee);
 867         }
 868         else
 869         {
 870             edat->s[i].ee = -1;
 871         }
 872     }
 873     sfree(eee);
 874 }
 875
 876 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
 877 {
 878     enerdata_t *esum;
 879     enerdat_t  *s;
 880     int         f, i;
 881     double      sum;
 882
 883     snew(esum, 1);
 884     *esum = *edat;
 885     snew(esum->s, 1);
 886     s = &esum->s[0];
 887     snew(s->ener, esum->nframes);
 888     snew(s->es, esum->nframes);
 889
 890     s->bExactStat = TRUE;
 891     s->slope      = 0;
 892     for (i = 0; i < nset; i++)
 893     {
 894         if (!edat->s[i].bExactStat)
 895         {
 896             s->bExactStat = FALSE;
 897         }
 898         s->slope += edat->s[i].slope;
 899     }
 900
 901     for (f = 0; f < edat->nframes; f++)
 902     {
 903         sum = 0;
 904         for (i = 0; i < nset; i++)
 905         {
 906             sum += edat->s[i].ener[f];
 907         }
 908         s->ener[f] = sum;
 909         sum        = 0;
 910         for (i = 0; i < nset; i++)
 911         {
 912             sum += edat->s[i].es[f].sum;
 913         }
 914         s->es[f].sum  = sum;
 915         s->es[f].sum2 = 0;
 916     }
 917
 918     calc_averages(1, esum, nbmin, nbmax);
 919
 920     return esum;
 921 }
 922
 923 static char *ee_pr(double ee, char *buf)
 924 {
 925     char   tmp[100];
 926     double rnd;
 927
 928     if (ee < 0)
 929     {
 930         sprintf(buf, "%s", "--");
 931     }
 932     else
 933     {
 934         /* Round to two decimals by printing. */
 935         sprintf(tmp, "%.1e", ee);
 936         sscanf(tmp, "%lf", &rnd);
 937         sprintf(buf, "%g", rnd);
 938     }
 939
 940     return buf;
 941 }
 942
 943 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
 944 {
 945 /* Remove the drift by performing a fit to y = ax+b.
 946    Uses 5 iterations. */
 947     int    i, j, k;
 948     double delta, d, sd, sd2;
 949
 950     edat->npoints = edat->nframes;
 951     edat->nsteps  = edat->nframes;
 952
 953     for (k = 0; (k < 5); k++)
 954     {
 955         for (i = 0; (i < nset); i++)
 956         {
 957             delta = edat->s[i].slope*dt;
 958
 959             if (NULL != debug)
 960             {
 961                 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
 962             }
 963
 964             for (j = 0; (j < edat->nframes); j++)
 965             {
 966                 edat->s[i].ener[j]   -= j*delta;
 967                 edat->s[i].es[j].sum  = 0;
 968                 edat->s[i].es[j].sum2 = 0;
 969             }
 970         }
 971         calc_averages(nset, edat, nbmin, nbmax);
 972     }
 973 }
 974
 975 static void calc_fluctuation_props(FILE *fp,
 976                                    gmx_bool bDriftCorr, real dt,
 977                                    int nset, int nmol,
 978                                    char **leg, enerdata_t *edat,
 979                                    int nbmin, int nbmax)
 980 {
 981     int    i, j;
 982     double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
 983     double NANO3;
 984     enum {
 985         eVol, eEnth, eTemp, eEtot, eNR
 986     };
 987     const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
 988     int         ii[eNR];
 989
 990     NANO3 = NANO*NANO*NANO;
 991     if (!bDriftCorr)
 992     {
 993         fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
 994                 "for spurious drift in the graphs. Note that this is not\n"
 995                 "a substitute for proper equilibration and sampling!\n");
 996     }
 997     else
 998     {
 999         remove_drift(nset, nbmin, nbmax, dt, edat);
1000     }
1001     for (i = 0; (i < eNR); i++)
1002     {
1003         for (ii[i] = 0; (ii[i] < nset &&
1004                          (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
1005         {
1006             ;
1007         }
1008 /*        if (ii[i] < nset)
1009             fprintf(fp,"Found %s data.\n",my_ener[i]);
1010  */ }
1011     /* Compute it all! */
1012     alpha = kappa = cp = dcp = cv = NOTSET;
1013
1014     /* Temperature */
1015     tt = NOTSET;
1016     if (ii[eTemp] < nset)
1017     {
1018         tt    = edat->s[ii[eTemp]].av;
1019     }
1020     /* Volume */
1021     vv = varv = NOTSET;
1022     if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1023     {
1024         vv    = edat->s[ii[eVol]].av*NANO3;
1025         varv  = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1026         kappa = (varv/vv)/(BOLTZMANN*tt);
1027     }
1028     /* Enthalpy */
1029     hh = varh = NOTSET;
1030     if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1031     {
1032         hh    = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1033         varh  = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1034         cp    = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1035     }
1036     /* Total energy */
1037     et = varet = NOTSET;
1038     if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1039     {
1040         /* Only compute cv in constant volume runs, which we can test
1041            by checking whether the enthalpy was computed.
1042          */
1043         et    = edat->s[ii[eEtot]].av;
1044         varet = sqr(edat->s[ii[eEtot]].rmsd);
1045         cv    = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1046     }
1047     /* Alpha, dcp */
1048     if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1049     {
1050         double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
1051         vh_sum = v_sum = h_sum = 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             v_sum  += v;
1057             h_sum  += h;
1058             vh_sum += (v*h);
1059         }
1060         vh_aver = vh_sum / edat->nframes;
1061         v_aver  = v_sum  / edat->nframes;
1062         h_aver  = h_sum  / edat->nframes;
1063         alpha   = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
1064         dcp     = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1065     }
1066
1067     if (tt != NOTSET)
1068     {
1069         if (nmol < 2)
1070         {
1071             fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1072                     nmol);
1073         }
1074         fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1075         fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1076         fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1077         fprintf(fp, "please use the g_dos program.\n\n");
1078         fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1079                 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1080
1081         if (debug != NULL)
1082         {
1083             if (varv != NOTSET)
1084             {
1085                 fprintf(fp, "varv  =  %10g (m^6)\n", varv*AVOGADRO/nmol);
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,
1135                          gmx_bool bVisco, const char *visfn, int nmol,
1136                          gmx_int64_t start_step, double start_t,
1137                          gmx_int64_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_int64_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);
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);
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             gmx_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         gmx_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_int64_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_int64_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_int64) ||
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         "[THISMODULE] extracts energy components or distance restraint",
1841         "data from an energy file. The user is prompted to interactively",
1842         "select the desired energy terms.[PAR]",
1843
1844         "Average, RMSD, and drift are calculated with full precision from the",
1845         "simulation (see printed manual). Drift is calculated by performing",
1846         "a least-squares fit of the data to a straight line. The reported total drift",
1847         "is the difference of the fit at the first and last point.",
1848         "An error estimate of the average is given based on a block averages",
1849         "over 5 blocks using the full-precision averages. The error estimate",
1850         "can be performed over multiple block lengths with the options",
1851         "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1852         "[BB]Note[bb] that in most cases the energy files contains averages over all",
1853         "MD steps, or over many more points than the number of frames in",
1854         "energy file. This makes the [THISMODULE] statistics output more accurate",
1855         "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1856         "file, the statistics mentioned above are simply over the single, per-frame",
1857         "energy values.[PAR]",
1858
1859         "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1860
1861         "Some fluctuation-dependent properties can be calculated provided",
1862         "the correct energy terms are selected, and that the command line option",
1863         "[TT]-fluct_props[tt] is given. The following properties",
1864         "will be computed:[BR]",
1865         "Property                        Energy terms needed[BR]",
1866         "---------------------------------------------------[BR]",
1867         "Heat capacity C[SUB]p[sub] (NPT sims):    Enthalpy, Temp     [BR]",
1868         "Heat capacity C[SUB]v[sub] (NVT sims):    Etot, Temp         [BR]",
1869         "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1870         "Isothermal compressibility:     Vol, Temp          [BR]",
1871         "Adiabatic bulk modulus:         Vol, Temp          [BR]",
1872         "---------------------------------------------------[BR]",
1873         "You always need to set the number of molecules [TT]-nmol[tt].",
1874         "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1875         "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]"
1876         "When the [TT]-viol[tt] option is set, the time averaged",
1877         "violations are plotted and the running time-averaged and",
1878         "instantaneous sum of violations are recalculated. Additionally",
1879         "running time-averaged and instantaneous distances between",
1880         "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1881
1882         "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1883         "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1884         "The first two options plot the orientation, the last three the",
1885         "deviations of the orientations from the experimental values.",
1886         "The options that end on an 'a' plot the average over time",
1887         "as a function of restraint. The options that end on a 't'",
1888         "prompt the user for restraint label numbers and plot the data",
1889         "as a function of time. Option [TT]-odr[tt] plots the RMS",
1890         "deviation as a function of restraint.",
1891         "When the run used time or ensemble averaged orientation restraints,",
1892         "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1893         "not ensemble-averaged orientations and deviations instead of",
1894         "the time and ensemble averages.[PAR]",
1895
1896         "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1897         "tensor for each orientation restraint experiment. With option",
1898         "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1899
1900         "Option [TT]-odh[tt] extracts and plots the free energy data",
1901         "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1902         "from the [TT]ener.edr[tt] file.[PAR]",
1903
1904         "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1905         "difference with an ideal gas state: [BR]",
1906         "  [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]",
1907         "  [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]",
1908         "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1909         "the average is over the ensemble (or time in a trajectory).",
1910         "Note that this is in principle",
1911         "only correct when averaging over the whole (Boltzmann) ensemble",
1912         "and using the potential energy. This also allows for an entropy",
1913         "estimate using:[BR]",
1914         "  [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]",
1915         "  [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",
1916         "[PAR]",
1917
1918         "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1919         "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1920         "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1921         "files, and the average is over the ensemble A. The running average",
1922         "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1923         "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1924
1925     };
1926     static gmx_bool    bSum    = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1927     static gmx_bool    bDp     = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1928     static int         skip    = 0, nmol = 1, nbmin = 5, nbmax = 5;
1929     static real        reftemp = 300.0, ezero = 0;
1930     t_pargs            pa[]    = {
1931         { "-fee",   FALSE, etBOOL,  {&bFee},
1932           "Do a free energy estimate" },
1933         { "-fetemp", FALSE, etREAL, {&reftemp},
1934           "Reference temperature for free energy calculation" },
1935         { "-zero", FALSE, etREAL, {&ezero},
1936           "Subtract a zero-point energy" },
1937         { "-sum",  FALSE, etBOOL, {&bSum},
1938           "Sum the energy terms selected rather than display them all" },
1939         { "-dp",   FALSE, etBOOL, {&bDp},
1940           "Print energies in high precision" },
1941         { "-nbmin", FALSE, etINT, {&nbmin},
1942           "Minimum number of blocks for error estimate" },
1943         { "-nbmax", FALSE, etINT, {&nbmax},
1944           "Maximum number of blocks for error estimate" },
1945         { "-mutot", FALSE, etBOOL, {&bMutot},
1946           "Compute the total dipole moment from the components" },
1947         { "-skip", FALSE, etINT,  {&skip},
1948           "Skip number of frames between data points" },
1949         { "-aver", FALSE, etBOOL, {&bPrAll},
1950           "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1951         { "-nmol", FALSE, etINT,  {&nmol},
1952           "Number of molecules in your sample: the energies are divided by this number" },
1953         { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1954           "Compute properties based on energy fluctuations, like heat capacity" },
1955         { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1956           "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1957         { "-fluc", FALSE, etBOOL, {&bFluct},
1958           "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1959         { "-orinst", FALSE, etBOOL, {&bOrinst},
1960           "Analyse instantaneous orientation data" },
1961         { "-ovec", FALSE, etBOOL, {&bOvec},
1962           "Also plot the eigenvectors with [TT]-oten[tt]" }
1963     };
1964     const char       * drleg[] = {
1965         "Running average",
1966         "Instantaneous"
1967     };
1968     static const char *setnm[] = {
1969         "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1970         "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1971         "Volume",  "Pressure"
1972     };
1973
1974     FILE              *out     = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1975     FILE              *fp_dhdl = NULL;
1976     FILE             **drout;
1977     ener_file_t        fp;
1978     int                timecheck = 0;
1979     gmx_mtop_t         mtop;
1980     gmx_localtop_t    *top = NULL;
1981     t_inputrec         ir;
1982     t_energy         **ee;
1983     enerdata_t         edat;
1984     gmx_enxnm_t       *enm = NULL;
1985     t_enxframe        *frame, *fr = NULL;
1986     int                cur = 0;
1987 #define NEXT (1-cur)
1988     int                nre, teller, teller_disre, nfr;
1989     gmx_int64_t        start_step;
1990     int                nor = 0, nex = 0, norfr = 0, enx_i = 0;
1991     real               start_t;
1992     real              *bounds  = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
1993     int               *index   = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
1994     int                nbounds = 0, npairs;
1995     gmx_bool           bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
1996     gmx_bool           bFoundStart, bCont, bEDR, bVisco;
1997     double             sum, sumaver, sumt, ener, dbl;
1998     double            *time = NULL;
1999     real               Vaver;
2000     int               *set     = NULL, i, j, k, nset, sss;
2001     gmx_bool          *bIsEner = NULL;
2002     char             **pairleg, **odtleg, **otenleg;
2003     char             **leg = NULL;
2004     char             **nms;
2005     char              *anm_j, *anm_k, *resnm_j, *resnm_k;
2006     int                resnr_j, resnr_k;
2007     const char        *orinst_sub = "@ subtitle \"instantaneous\"\n";
2008     char               buf[256];
2009     output_env_t       oenv;
2010     t_enxblock        *blk       = NULL;
2011     t_enxblock        *blk_disre = NULL;
2012     int                ndisre    = 0;
2013     int                dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2014
2015     t_filenm           fnm[] = {
2016         { efEDR, "-f",    NULL,      ffREAD  },
2017         { efEDR, "-f2",   NULL,      ffOPTRD },
2018         { efTPX, "-s",    NULL,      ffOPTRD },
2019         { efXVG, "-o",    "energy",  ffWRITE },
2020         { efXVG, "-viol", "violaver", ffOPTWR },
2021         { efXVG, "-pairs", "pairs",   ffOPTWR },
2022         { efXVG, "-ora",  "orienta", ffOPTWR },
2023         { efXVG, "-ort",  "orientt", ffOPTWR },
2024         { efXVG, "-oda",  "orideva", ffOPTWR },
2025         { efXVG, "-odr",  "oridevr", ffOPTWR },
2026         { efXVG, "-odt",  "oridevt", ffOPTWR },
2027         { efXVG, "-oten", "oriten",  ffOPTWR },
2028         { efXVG, "-corr", "enecorr", ffOPTWR },
2029         { efXVG, "-vis",  "visco",   ffOPTWR },
2030         { efXVG, "-ravg", "runavgdf", ffOPTWR },
2031         { efXVG, "-odh",  "dhdl", ffOPTWR }
2032     };
2033 #define NFILE asize(fnm)
2034     int                npargs;
2035     t_pargs           *ppa;
2036
2037     npargs = asize(pa);
2038     ppa    = add_acf_pargs(&npargs, pa);
2039     if (!parse_common_args(&argc, argv,
2040                            PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
2041                            NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
2042     {
2043         return 0;
2044     }
2045
2046     bDRAll = opt2bSet("-pairs", NFILE, fnm);
2047     bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2048     bORA   = opt2bSet("-ora", NFILE, fnm);
2049     bORT   = opt2bSet("-ort", NFILE, fnm);
2050     bODA   = opt2bSet("-oda", NFILE, fnm);
2051     bODR   = opt2bSet("-odr", NFILE, fnm);
2052     bODT   = opt2bSet("-odt", NFILE, fnm);
2053     bORIRE = bORA || bORT || bODA || bODR || bODT;
2054     bOTEN  = opt2bSet("-oten", NFILE, fnm);
2055     bDHDL  = opt2bSet("-odh", NFILE, fnm);
2056
2057     nset = 0;
2058
2059     snew(frame, 2);
2060     fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2061     do_enxnms(fp, &nre, &enm);
2062
2063     Vaver = -1;
2064
2065     bVisco = opt2bSet("-vis", NFILE, fnm);
2066
2067     if ((!bDisRe) && (!bDHDL))
2068     {
2069         if (bVisco)
2070         {
2071             nset = asize(setnm);
2072             snew(set, nset);
2073             /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2074             for (j = 0; j < nset; j++)
2075             {
2076                 for (i = 0; i < nre; i++)
2077                 {
2078                     if (strstr(enm[i].name, setnm[j]))
2079                     {
2080                         set[j] = i;
2081                         break;
2082                     }
2083                 }
2084                 if (i == nre)
2085                 {
2086                     if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2087                     {
2088                         printf("Enter the box volume (" unit_volume "): ");
2089                         if (1 != scanf("%lf", &dbl))
2090                         {
2091                             gmx_fatal(FARGS, "Error reading user input");
2092                         }
2093                         Vaver = dbl;
2094                     }
2095                     else
2096                     {
2097                         gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2098                                   setnm[j]);
2099                     }
2100                 }
2101             }
2102         }
2103         else
2104         {
2105             set = select_by_name(nre, enm, &nset);
2106         }
2107         /* Print all the different units once */
2108         sprintf(buf, "(%s)", enm[set[0]].unit);
2109         for (i = 1; i < nset; i++)
2110         {
2111             for (j = 0; j < i; j++)
2112             {
2113                 if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2114                 {
2115                     break;
2116                 }
2117             }
2118             if (j == i)
2119             {
2120                 strcat(buf, ", (");
2121                 strcat(buf, enm[set[i]].unit);
2122                 strcat(buf, ")");
2123             }
2124         }
2125         out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf,
2126                        oenv);
2127
2128         snew(leg, nset+1);
2129         for (i = 0; (i < nset); i++)
2130         {
2131             leg[i] = enm[set[i]].name;
2132         }
2133         if (bSum)
2134         {
2135             leg[nset] = strdup("Sum");
2136             xvgr_legend(out, nset+1, (const char**)leg, oenv);
2137         }
2138         else
2139         {
2140             xvgr_legend(out, nset, (const char**)leg, oenv);
2141         }
2142
2143         snew(bIsEner, nset);
2144         for (i = 0; (i < nset); i++)
2145         {
2146             bIsEner[i] = FALSE;
2147             for (j = 0; (j <= F_ETOT); j++)
2148             {
2149                 bIsEner[i] = bIsEner[i] ||
2150                     (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2151             }
2152         }
2153
2154         if (bPrAll && nset > 1)
2155         {
2156             gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2157         }
2158
2159         time = NULL;
2160
2161         if (bORIRE || bOTEN)
2162         {
2163             get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2164         }
2165
2166         if (bORIRE)
2167         {
2168             if (bOrinst)
2169             {
2170                 enx_i = enxORI;
2171             }
2172             else
2173             {
2174                 enx_i = enxOR;
2175             }
2176
2177             if (bORA || bODA)
2178             {
2179                 snew(orient, nor);
2180             }
2181             if (bODR)
2182             {
2183                 snew(odrms, nor);
2184             }
2185             if (bORT || bODT)
2186             {
2187                 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2188                 fprintf(stderr, "End your selection with 0\n");
2189                 j     = -1;
2190                 orsel = NULL;
2191                 do
2192                 {
2193                     j++;
2194                     srenew(orsel, j+1);
2195                     if (1 != scanf("%d", &(orsel[j])))
2196                     {
2197                         gmx_fatal(FARGS, "Error reading user input");
2198                     }
2199                 }
2200                 while (orsel[j] > 0);
2201                 if (orsel[0] == -1)
2202                 {
2203                     fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2204                     norsel = nor;
2205                     srenew(orsel, nor);
2206                     for (i = 0; i < nor; i++)
2207                     {
2208                         orsel[i] = i;
2209                     }
2210                 }
2211                 else
2212                 {
2213                     /* Build the selection */
2214                     norsel = 0;
2215                     for (i = 0; i < j; i++)
2216                     {
2217                         for (k = 0; k < nor; k++)
2218                         {
2219                             if (or_label[k] == orsel[i])
2220                             {
2221                                 orsel[norsel] = k;
2222                                 norsel++;
2223                                 break;
2224                             }
2225                         }
2226                         if (k == nor)
2227                         {
2228                             fprintf(stderr, "Orientation restraint label %d not found\n",
2229                                     orsel[i]);
2230                         }
2231                     }
2232                 }
2233                 snew(odtleg, norsel);
2234                 for (i = 0; i < norsel; i++)
2235                 {
2236                     snew(odtleg[i], 256);
2237                     sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2238                 }
2239                 if (bORT)
2240                 {
2241                     fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2242                                     "Time (ps)", "", oenv);
2243                     if (bOrinst)
2244                     {
2245                         fprintf(fort, "%s", orinst_sub);
2246                     }
2247                     xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2248                 }
2249                 if (bODT)
2250                 {
2251                     fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2252                                     "Orientation restraint deviation",
2253                                     "Time (ps)", "", oenv);
2254                     if (bOrinst)
2255                     {
2256                         fprintf(fodt, "%s", orinst_sub);
2257                     }
2258                     xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2259                 }
2260             }
2261         }
2262         if (bOTEN)
2263         {
2264             foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2265                              "Order tensor", "Time (ps)", "", oenv);
2266             snew(otenleg, bOvec ? nex*12 : nex*3);
2267             for (i = 0; i < nex; i++)
2268             {
2269                 for (j = 0; j < 3; j++)
2270                 {
2271                     sprintf(buf, "eig%d", j+1);
2272                     otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
2273                 }
2274                 if (bOvec)
2275                 {
2276                     for (j = 0; j < 9; j++)
2277                     {
2278                         sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2279                         otenleg[12*i+3+j] = strdup(buf);
2280                     }
2281                 }
2282             }
2283             xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2284         }
2285     }
2286     else if (bDisRe)
2287     {
2288         nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
2289                              &mtop, &top, &ir);
2290         snew(violaver, npairs);
2291         out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2292                        "Time (ps)", "nm", oenv);
2293         xvgr_legend(out, 2, drleg, oenv);
2294         if (bDRAll)
2295         {
2296             fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2297                                 "Time (ps)", "Distance (nm)", oenv);
2298             if (output_env_get_print_xvgr_codes(oenv))
2299             {
2300                 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2301                         ir.dr_tau);
2302             }
2303         }
2304     }
2305     else if (bDHDL)
2306     {
2307         get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
2308     }
2309
2310     /* Initiate energies and set them to zero */
2311     edat.nsteps  = 0;
2312     edat.npoints = 0;
2313     edat.nframes = 0;
2314     edat.step    = NULL;
2315     edat.steps   = NULL;
2316     edat.points  = NULL;
2317     snew(edat.s, nset);
2318
2319     /* Initiate counters */
2320     teller       = 0;
2321     teller_disre = 0;
2322     bFoundStart  = FALSE;
2323     start_step   = 0;
2324     start_t      = 0;
2325     do
2326     {
2327         /* This loop searches for the first frame (when -b option is given),
2328          * or when this has been found it reads just one energy frame
2329          */
2330         do
2331         {
2332             bCont = do_enx(fp, &(frame[NEXT]));
2333             if (bCont)
2334             {
2335                 timecheck = check_times(frame[NEXT].t);
2336             }
2337         }
2338         while (bCont && (timecheck < 0));
2339
2340         if ((timecheck == 0) && bCont)
2341         {
2342             /* We read a valid frame, so we can use it */
2343             fr = &(frame[NEXT]);
2344
2345             if (fr->nre > 0)
2346             {
2347                 /* The frame contains energies, so update cur */
2348                 cur  = NEXT;
2349
2350                 if (edat.nframes % 1000 == 0)
2351                 {
2352                     srenew(edat.step, edat.nframes+1000);
2353                     memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2354                     srenew(edat.steps, edat.nframes+1000);
2355                     memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2356                     srenew(edat.points, edat.nframes+1000);
2357                     memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2358
2359                     for (i = 0; i < nset; i++)
2360                     {
2361                         srenew(edat.s[i].ener, edat.nframes+1000);
2362                         memset(&(edat.s[i].ener[edat.nframes]), 0,
2363                                1000*sizeof(edat.s[i].ener[0]));
2364                         srenew(edat.s[i].es, edat.nframes+1000);
2365                         memset(&(edat.s[i].es[edat.nframes]), 0,
2366                                1000*sizeof(edat.s[i].es[0]));
2367                     }
2368                 }
2369
2370                 nfr            = edat.nframes;
2371                 edat.step[nfr] = fr->step;
2372
2373                 if (!bFoundStart)
2374                 {
2375                     bFoundStart = TRUE;
2376                     /* Initiate the previous step data */
2377                     start_step = fr->step;
2378                     start_t    = fr->t;
2379                     /* Initiate the energy sums */
2380                     edat.steps[nfr]  = 1;
2381                     edat.points[nfr] = 1;
2382                     for (i = 0; i < nset; i++)
2383                     {
2384                         sss                    = set[i];
2385                         edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2386                         edat.s[i].es[nfr].sum2 = 0;
2387                     }
2388                     edat.nsteps  = 1;
2389                     edat.npoints = 1;
2390                 }
2391                 else
2392                 {
2393                     edat.steps[nfr] = fr->nsteps;
2394                     {
2395                         if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2396                         {
2397                             if (fr->nsum <= 1)
2398                             {
2399                                 edat.points[nfr] = 1;
2400                                 for (i = 0; i < nset; i++)
2401                                 {
2402                                     sss                    = set[i];
2403                                     edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2404                                     edat.s[i].es[nfr].sum2 = 0;
2405                                 }
2406                                 edat.npoints += 1;
2407                             }
2408                             else
2409                             {
2410                                 edat.points[nfr] = fr->nsum;
2411                                 for (i = 0; i < nset; i++)
2412                                 {
2413                                     sss                    = set[i];
2414                                     edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2415                                     edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2416                                 }
2417                                 edat.npoints += fr->nsum;
2418                             }
2419                         }
2420                         else
2421                         {
2422                             /* The interval does not match fr->nsteps:
2423                              * can not do exact averages.
2424                              */
2425                             edat.npoints = 0;
2426                         }
2427                         edat.nsteps = fr->step - start_step + 1;
2428                     }
2429                 }
2430                 for (i = 0; i < nset; i++)
2431                 {
2432                     edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2433                 }
2434             }
2435             /*
2436              * Define distance restraint legends. Can only be done after
2437              * the first frame has been read... (Then we know how many there are)
2438              */
2439             blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2440             if (bDisRe && bDRAll && !leg && blk_disre)
2441             {
2442                 t_iatom   *fa;
2443                 t_iparams *ip;
2444
2445                 fa = top->idef.il[F_DISRES].iatoms;
2446                 ip = top->idef.iparams;
2447                 if (blk_disre->nsub != 2 ||
2448                     (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2449                 {
2450                     gmx_incons("Number of disre sub-blocks not equal to 2");
2451                 }
2452
2453                 ndisre = blk_disre->sub[0].nr;
2454                 if (ndisre != top->idef.il[F_DISRES].nr/3)
2455                 {
2456                     gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2457                               ndisre, top->idef.il[F_DISRES].nr/3);
2458                 }
2459                 snew(pairleg, ndisre);
2460                 for (i = 0; i < ndisre; i++)
2461                 {
2462                     snew(pairleg[i], 30);
2463                     j = fa[3*i+1];
2464                     k = fa[3*i+2];
2465                     gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2466                     gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2467                     sprintf(pairleg[i], "%d %s %d %s (%d)",
2468                             resnr_j, anm_j, resnr_k, anm_k,
2469                             ip[fa[3*i]].disres.label);
2470                 }
2471                 set = select_it(ndisre, pairleg, &nset);
2472                 snew(leg, 2*nset);
2473                 for (i = 0; (i < nset); i++)
2474                 {
2475                     snew(leg[2*i], 32);
2476                     sprintf(leg[2*i],  "a %s", pairleg[set[i]]);
2477                     snew(leg[2*i+1], 32);
2478                     sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2479                 }
2480                 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2481             }
2482
2483             /*
2484              * Store energies for analysis afterwards...
2485              */
2486             if (!bDisRe && !bDHDL && (fr->nre > 0))
2487             {
2488                 if (edat.nframes % 1000 == 0)
2489                 {
2490                     srenew(time, edat.nframes+1000);
2491                 }
2492                 time[edat.nframes] = fr->t;
2493                 edat.nframes++;
2494             }
2495             /*
2496              * Printing time, only when we do not want to skip frames
2497              */
2498             if (!skip || teller % skip == 0)
2499             {
2500                 if (bDisRe)
2501                 {
2502                     /*******************************************
2503                      * D I S T A N C E   R E S T R A I N T S
2504                      *******************************************/
2505                     if (ndisre > 0)
2506                     {
2507  #ifndef GMX_DOUBLE
2508                         float  *disre_rt     =     blk_disre->sub[0].fval;
2509                         float  *disre_rm3tav = blk_disre->sub[1].fval;
2510  #else
2511                         double *disre_rt     =     blk_disre->sub[0].dval;
2512                         double *disre_rm3tav = blk_disre->sub[1].dval;
2513  #endif
2514
2515                         print_time(out, fr->t);
2516                         if (violaver == NULL)
2517                         {
2518                             snew(violaver, ndisre);
2519                         }
2520
2521                         /* Subtract bounds from distances, to calculate violations */
2522                         calc_violations(disre_rt, disre_rm3tav,
2523                                         nbounds, pair, bounds, violaver, &sumt, &sumaver);
2524
2525                         fprintf(out, "  %8.4f  %8.4f\n", sumaver, sumt);
2526                         if (bDRAll)
2527                         {
2528                             print_time(fp_pairs, fr->t);
2529                             for (i = 0; (i < nset); i++)
2530                             {
2531                                 sss = set[i];
2532                                 fprintf(fp_pairs, "  %8.4f", mypow(disre_rm3tav[sss], minthird));
2533                                 fprintf(fp_pairs, "  %8.4f", disre_rt[sss]);
2534                             }
2535                             fprintf(fp_pairs, "\n");
2536                         }
2537                         teller_disre++;
2538                     }
2539                 }
2540                 else if (bDHDL)
2541                 {
2542                     do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2543                 }
2544
2545                 /*******************************************
2546                  * E N E R G I E S
2547                  *******************************************/
2548                 else
2549                 {
2550                     if (fr->nre > 0)
2551                     {
2552                         if (bPrAll)
2553                         {
2554                             /* We skip frames with single points (usually only the first frame),
2555                              * since they would result in an average plot with outliers.
2556                              */
2557                             if (fr->nsum > 1)
2558                             {
2559                                 print_time(out, fr->t);
2560                                 print1(out, bDp, fr->ener[set[0]].e);
2561                                 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2562                                 print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
2563                                 fprintf(out, "\n");
2564                             }
2565                         }
2566                         else
2567                         {
2568                             print_time(out, fr->t);
2569                             if (bSum)
2570                             {
2571                                 sum = 0;
2572                                 for (i = 0; i < nset; i++)
2573                                 {
2574                                     sum += fr->ener[set[i]].e;
2575                                 }
2576                                 print1(out, bDp, sum/nmol-ezero);
2577                             }
2578                             else
2579                             {
2580                                 for (i = 0; (i < nset); i++)
2581                                 {
2582                                     if (bIsEner[i])
2583                                     {
2584                                         print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2585                                     }
2586                                     else
2587                                     {
2588                                         print1(out, bDp, fr->ener[set[i]].e);
2589                                     }
2590                                 }
2591                             }
2592                             fprintf(out, "\n");
2593                         }
2594                     }
2595                     blk = find_block_id_enxframe(fr, enx_i, NULL);
2596                     if (bORIRE && blk)
2597                     {
2598 #ifndef GMX_DOUBLE
2599                         xdr_datatype dt = xdr_datatype_float;
2600 #else
2601                         xdr_datatype dt = xdr_datatype_double;
2602 #endif
2603                         real        *vals;
2604
2605                         if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2606                         {
2607                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2608                         }
2609 #ifndef GMX_DOUBLE
2610                         vals = blk->sub[0].fval;
2611 #else
2612                         vals = blk->sub[0].dval;
2613 #endif
2614
2615                         if (blk->sub[0].nr != (size_t)nor)
2616                         {
2617                             gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2618                         }
2619                         if (bORA || bODA)
2620                         {
2621                             for (i = 0; i < nor; i++)
2622                             {
2623                                 orient[i] += vals[i];
2624                             }
2625                         }
2626                         if (bODR)
2627                         {
2628                             for (i = 0; i < nor; i++)
2629                             {
2630                                 odrms[i] += sqr(vals[i]-oobs[i]);
2631                             }
2632                         }
2633                         if (bORT)
2634                         {
2635                             fprintf(fort, "  %10f", fr->t);
2636                             for (i = 0; i < norsel; i++)
2637                             {
2638                                 fprintf(fort, " %g", vals[orsel[i]]);
2639                             }
2640                             fprintf(fort, "\n");
2641                         }
2642                         if (bODT)
2643                         {
2644                             fprintf(fodt, "  %10f", fr->t);
2645                             for (i = 0; i < norsel; i++)
2646                             {
2647                                 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2648                             }
2649                             fprintf(fodt, "\n");
2650                         }
2651                         norfr++;
2652                     }
2653                     blk = find_block_id_enxframe(fr, enxORT, NULL);
2654                     if (bOTEN && blk)
2655                     {
2656 #ifndef GMX_DOUBLE
2657                         xdr_datatype dt = xdr_datatype_float;
2658 #else
2659                         xdr_datatype dt = xdr_datatype_double;
2660 #endif
2661                         real        *vals;
2662
2663                         if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2664                         {
2665                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2666                         }
2667 #ifndef GMX_DOUBLE
2668                         vals = blk->sub[0].fval;
2669 #else
2670                         vals = blk->sub[0].dval;
2671 #endif
2672
2673                         if (blk->sub[0].nr != (size_t)(nex*12))
2674                         {
2675                             gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2676                                       blk->sub[0].nr/12, nex);
2677                         }
2678                         fprintf(foten, "  %10f", fr->t);
2679                         for (i = 0; i < nex; i++)
2680                         {
2681                             for (j = 0; j < (bOvec ? 12 : 3); j++)
2682                             {
2683                                 fprintf(foten, " %g", vals[i*12+j]);
2684                             }
2685                         }
2686                         fprintf(foten, "\n");
2687                     }
2688                 }
2689             }
2690             teller++;
2691         }
2692     }
2693     while (bCont && (timecheck == 0));
2694
2695     fprintf(stderr, "\n");
2696     close_enx(fp);
2697     if (out)
2698     {
2699         gmx_ffclose(out);
2700     }
2701
2702     if (bDRAll)
2703     {
2704         gmx_ffclose(fp_pairs);
2705     }
2706
2707     if (bORT)
2708     {
2709         gmx_ffclose(fort);
2710     }
2711     if (bODT)
2712     {
2713         gmx_ffclose(fodt);
2714     }
2715     if (bORA)
2716     {
2717         out = xvgropen(opt2fn("-ora", NFILE, fnm),
2718                        "Average calculated orientations",
2719                        "Restraint label", "", oenv);
2720         if (bOrinst)
2721         {
2722             fprintf(out, "%s", orinst_sub);
2723         }
2724         for (i = 0; i < nor; i++)
2725         {
2726             fprintf(out, "%5d  %g\n", or_label[i], orient[i]/norfr);
2727         }
2728         gmx_ffclose(out);
2729     }
2730     if (bODA)
2731     {
2732         out = xvgropen(opt2fn("-oda", NFILE, fnm),
2733                        "Average restraint deviation",
2734                        "Restraint label", "", oenv);
2735         if (bOrinst)
2736         {
2737             fprintf(out, "%s", orinst_sub);
2738         }
2739         for (i = 0; i < nor; i++)
2740         {
2741             fprintf(out, "%5d  %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2742         }
2743         gmx_ffclose(out);
2744     }
2745     if (bODR)
2746     {
2747         out = xvgropen(opt2fn("-odr", NFILE, fnm),
2748                        "RMS orientation restraint deviations",
2749                        "Restraint label", "", oenv);
2750         if (bOrinst)
2751         {
2752             fprintf(out, "%s", orinst_sub);
2753         }
2754         for (i = 0; i < nor; i++)
2755         {
2756             fprintf(out, "%5d  %g\n", or_label[i], sqrt(odrms[i]/norfr));
2757         }
2758         gmx_ffclose(out);
2759     }
2760     if (bOTEN)
2761     {
2762         gmx_ffclose(foten);
2763     }
2764
2765     if (bDisRe)
2766     {
2767         analyse_disre(opt2fn("-viol", NFILE, fnm),
2768                       teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2769     }
2770     else if (bDHDL)
2771     {
2772         if (fp_dhdl)
2773         {
2774             gmx_ffclose(fp_dhdl);
2775             printf("\n\nWrote %d lambda values with %d samples as ",
2776                    dh_lambdas, dh_samples);
2777             if (dh_hists > 0)
2778             {
2779                 printf("%d dH histograms ", dh_hists);
2780             }
2781             if (dh_blocks > 0)
2782             {
2783                 printf("%d dH data blocks ", dh_blocks);
2784             }
2785             printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2786
2787         }
2788         else
2789         {
2790             gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2791         }
2792
2793     }
2794     else
2795     {
2796         double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2797         analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2798                      bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa),
2799                      bVisco, opt2fn("-vis", NFILE, fnm),
2800                      nmol,
2801                      start_step, start_t, frame[cur].step, frame[cur].t,
2802                      time, reftemp, &edat,
2803                      nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2804                      oenv);
2805         if (bFluctProps)
2806         {
2807             calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2808                                    nbmin, nbmax);
2809         }
2810     }
2811     if (opt2bSet("-f2", NFILE, fnm))
2812     {
2813         fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2814             reftemp, nset, set, leg, &edat, time, oenv);
2815     }
2816
2817     {
2818         const char *nxy = "-nxy";
2819
2820         do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2821         do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2822         do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2823         do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2824         do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2825         do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2826         do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2827         do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2828         do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);
2829     }
2830
2831     return 0;
2832 }