src/gromacs/gmxana/gmx_energy.c

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