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