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