src/tools/gmx_energy.c

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