src/tools/gmx_energy.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <string.h>
  43 #include <math.h>
  44 #include <ctype.h>
  45
  46 #include "typedefs.h"
  47 #include "gmx_fatal.h"
  48 #include "vec.h"
  49 #include "string2.h"
  50 #include "smalloc.h"
  51 #include "enxio.h"
  52 #include "statutil.h"
  53 #include "names.h"
  54 #include "copyrite.h"
  55 #include "macros.h"
  56 #include "xvgr.h"
  57 #include "gstat.h"
  58 #include "physics.h"
  59 #include "tpxio.h"
  60 #include "viewit.h"
  61 #include "mtop_util.h"
  62 #include "gmx_ana.h"
  63 #include "mdebin.h"
  64
  65 static real       minthird = -1.0/3.0, minsixth = -1.0/6.0;
  66
  67 typedef struct {
  68     real sum;
  69     real sum2;
  70 } exactsum_t;
  71
  72 typedef struct {
  73     real       *ener;
  74     exactsum_t *es;
  75     gmx_bool    bExactStat;
  76     double      av;
  77     double      rmsd;
  78     double      ee;
  79     double      slope;
  80 } enerdat_t;
  81
  82 typedef struct {
  83     gmx_large_int_t  nsteps;
  84     gmx_large_int_t  npoints;
  85     int              nframes;
  86     int             *step;
  87     int             *steps;
  88     int             *points;
  89     enerdat_t       *s;
  90 } enerdata_t;
  91
  92 static double mypow(double x, double y)
  93 {
  94     if (x > 0)
  95     {
  96         return pow(x, y);
  97     }
  98     else
  99     {
 100         return 0.0;
 101     }
 102 }
 103
 104 static int *select_it(int nre, char *nm[], int *nset)
 105 {
 106     gmx_bool *bE;
 107     int       n, k, j, i;
 108     int      *set;
 109     gmx_bool  bVerbose = TRUE;
 110
 111     if ((getenv("VERBOSE")) != NULL)
 112     {
 113         bVerbose = FALSE;
 114     }
 115
 116     fprintf(stderr, "Select the terms you want from the following list\n");
 117     fprintf(stderr, "End your selection with 0\n");
 118
 119     if (bVerbose)
 120     {
 121         for (k = 0; (k < nre); )
 122         {
 123             for (j = 0; (j < 4) && (k < nre); j++, k++)
 124             {
 125                 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
 126             }
 127             fprintf(stderr, "\n");
 128         }
 129     }
 130
 131     snew(bE, nre);
 132     do
 133     {
 134         if (1 != scanf("%d", &n))
 135         {
 136             gmx_fatal(FARGS, "Error reading user input");
 137         }
 138         if ((n > 0) && (n <= nre))
 139         {
 140             bE[n-1] = TRUE;
 141         }
 142     }
 143     while (n != 0);
 144
 145     snew(set, nre);
 146     for (i = (*nset) = 0; (i < nre); i++)
 147     {
 148         if (bE[i])
 149         {
 150             set[(*nset)++] = i;
 151         }
 152     }
 153
 154     sfree(bE);
 155
 156     return set;
 157 }
 158
 159 static void chomp(char *buf)
 160 {
 161     int len = strlen(buf);
 162
 163     while ((len > 0) && (buf[len-1] == '\n'))
 164     {
 165         buf[len-1] = '\0';
 166         len--;
 167     }
 168 }
 169
 170 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
 171 {
 172     gmx_bool   *bE;
 173     int         n, k, kk, j, i, nmatch, nind, nss;
 174     int        *set;
 175     gmx_bool    bEOF, bVerbose = TRUE, bLong = FALSE;
 176     char       *ptr, buf[STRLEN];
 177     const char *fm4   = "%3d  %-14s";
 178     const char *fm2   = "%3d  %-34s";
 179     char      **newnm = NULL;
 180
 181     if ((getenv("VERBOSE")) != NULL)
 182     {
 183         bVerbose = FALSE;
 184     }
 185
 186     fprintf(stderr, "\n");
 187     fprintf(stderr, "Select the terms you want from the following list by\n");
 188     fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
 189     fprintf(stderr, "End your selection with an empty line or a zero.\n");
 190     fprintf(stderr, "-------------------------------------------------------------------\n");
 191
 192     snew(newnm, nre);
 193     j = 0;
 194     for (k = 0; k < nre; k++)
 195     {
 196         newnm[k] = strdup(nm[k].name);
 197         /* Insert dashes in all the names */
 198         while ((ptr = strchr(newnm[k], ' ')) != NULL)
 199         {
 200             *ptr = '-';
 201         }
 202         if (bVerbose)
 203         {
 204             if (j == 0)
 205             {
 206                 if (k > 0)
 207                 {
 208                     fprintf(stderr, "\n");
 209                 }
 210                 bLong = FALSE;
 211                 for (kk = k; kk < k+4; kk++)
 212                 {
 213                     if (kk < nre && strlen(nm[kk].name) > 14)
 214                     {
 215                         bLong = TRUE;
 216                     }
 217                 }
 218             }
 219             else
 220             {
 221                 fprintf(stderr, " ");
 222             }
 223             if (!bLong)
 224             {
 225                 fprintf(stderr, fm4, k+1, newnm[k]);
 226                 j++;
 227                 if (j == 4)
 228                 {
 229                     j = 0;
 230                 }
 231             }
 232             else
 233             {
 234                 fprintf(stderr, fm2, k+1, newnm[k]);
 235                 j++;
 236                 if (j == 2)
 237                 {
 238                     j = 0;
 239                 }
 240             }
 241         }
 242     }
 243     if (bVerbose)
 244     {
 245         fprintf(stderr, "\n\n");
 246     }
 247
 248     snew(bE, nre);
 249
 250     bEOF = FALSE;
 251     while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
 252     {
 253         /* Remove newlines */
 254         chomp(buf);
 255
 256         /* Remove spaces */
 257         trim(buf);
 258
 259         /* Empty line means end of input */
 260         bEOF = (strlen(buf) == 0);
 261         if (!bEOF)
 262         {
 263             ptr = buf;
 264             do
 265             {
 266                 if (!bEOF)
 267                 {
 268                     /* First try to read an integer */
 269                     nss   = sscanf(ptr, "%d", &nind);
 270                     if (nss == 1)
 271                     {
 272                         /* Zero means end of input */
 273                         if (nind == 0)
 274                         {
 275                             bEOF = TRUE;
 276                         }
 277                         else if ((1 <= nind) && (nind <= nre))
 278                         {
 279                             bE[nind-1] = TRUE;
 280                         }
 281                         else
 282                         {
 283                             fprintf(stderr, "number %d is out of range\n", nind);
 284                         }
 285                     }
 286                     else
 287                     {
 288                         /* Now try to read a string */
 289                         i      = strlen(ptr);
 290                         nmatch = 0;
 291                         for (nind = 0; nind < nre; nind++)
 292                         {
 293                             if (gmx_strcasecmp(newnm[nind], ptr) == 0)
 294                             {
 295                                 bE[nind] = TRUE;
 296                                 nmatch++;
 297                             }
 298                         }
 299                         if (nmatch == 0)
 300                         {
 301                             i      = strlen(ptr);
 302                             nmatch = 0;
 303                             for (nind = 0; nind < nre; nind++)
 304                             {
 305                                 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
 306                                 {
 307                                     bE[nind] = TRUE;
 308                                     nmatch++;
 309                                 }
 310                             }
 311                             if (nmatch == 0)
 312                             {
 313                                 fprintf(stderr, "String '%s' does not match anything\n", ptr);
 314                             }
 315                         }
 316                     }
 317                 }
 318                 /* Look for the first space, and remove spaces from there */
 319                 if ((ptr = strchr(ptr, ' ')) != NULL)
 320                 {
 321                     trim(ptr);
 322                 }
 323             }
 324             while (!bEOF && (ptr && (strlen(ptr) > 0)));
 325         }
 326     }
 327
 328     snew(set, nre);
 329     for (i = (*nset) = 0; (i < nre); i++)
 330     {
 331         if (bE[i])
 332         {
 333             set[(*nset)++] = i;
 334         }
 335     }
 336
 337     sfree(bE);
 338
 339     if (*nset == 0)
 340     {
 341         gmx_fatal(FARGS, "No energy terms selected");
 342     }
 343
 344     for (i = 0; (i < nre); i++)
 345     {
 346         sfree(newnm[i]);
 347     }
 348     sfree(newnm);
 349
 350     return set;
 351 }
 352
 353 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
 354 {
 355     gmx_mtop_t  mtop;
 356     int         natoms;
 357     t_iatom    *iatom;
 358     matrix      box;
 359
 360     /* all we need is the ir to be able to write the label */
 361     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
 362 }
 363
 364 static void get_orires_parms(const char *topnm,
 365                              int *nor, int *nex, int **label, real **obs)
 366 {
 367     gmx_mtop_t      mtop;
 368     gmx_localtop_t *top;
 369     t_inputrec      ir;
 370     t_iparams      *ip;
 371     int             natoms, i;
 372     t_iatom        *iatom;
 373     int             nb;
 374     matrix          box;
 375
 376     read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
 377     top = gmx_mtop_generate_local_top(&mtop, &ir);
 378
 379     ip       = top->idef.iparams;
 380     iatom    = top->idef.il[F_ORIRES].iatoms;
 381
 382     /* Count how many distance restraint there are... */
 383     nb = top->idef.il[F_ORIRES].nr;
 384     if (nb == 0)
 385     {
 386         gmx_fatal(FARGS, "No orientation restraints in topology!\n");
 387     }
 388
 389     *nor = nb/3;
 390     *nex = 0;
 391     snew(*label, *nor);
 392     snew(*obs, *nor);
 393     for (i = 0; i < nb; i += 3)
 394     {
 395         (*label)[i/3] = ip[iatom[i]].orires.label;
 396         (*obs)[i/3]   = ip[iatom[i]].orires.obs;
 397         if (ip[iatom[i]].orires.ex >= *nex)
 398         {
 399             *nex = ip[iatom[i]].orires.ex+1;
 400         }
 401     }
 402     fprintf(stderr, "Found %d orientation restraints with %d experiments",
 403             *nor, *nex);
 404 }
 405
 406 static int get_bounds(const char *topnm,
 407                       real **bounds, int **index, int **dr_pair, int *npairs,
 408                       gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
 409 {
 410     gmx_localtop_t *top;
 411     t_functype     *functype;
 412     t_iparams      *ip;
 413     int             natoms, i, j, k, type, ftype, natom;
 414     t_ilist        *disres;
 415     t_iatom        *iatom;
 416     real           *b;
 417     int            *ind, *pair;
 418     int             nb, label1;
 419     matrix          box;
 420
 421     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
 422     snew(*ltop, 1);
 423     top   = gmx_mtop_generate_local_top(mtop, ir);
 424     *ltop = top;
 425
 426     functype = top->idef.functype;
 427     ip       = top->idef.iparams;
 428
 429     /* Count how many distance restraint there are... */
 430     nb = top->idef.il[F_DISRES].nr;
 431     if (nb == 0)
 432     {
 433         gmx_fatal(FARGS, "No distance restraints in topology!\n");
 434     }
 435
 436     /* Allocate memory */
 437     snew(b, nb);
 438     snew(ind, nb);
 439     snew(pair, nb+1);
 440
 441     /* Fill the bound array */
 442     nb = 0;
 443     for (i = 0; (i < top->idef.ntypes); i++)
 444     {
 445         ftype = functype[i];
 446         if (ftype == F_DISRES)
 447         {
 448
 449             label1  = ip[i].disres.label;
 450             b[nb]   = ip[i].disres.up1;
 451             ind[nb] = label1;
 452             nb++;
 453         }
 454     }
 455     *bounds = b;
 456
 457     /* Fill the index array */
 458     label1  = -1;
 459     disres  = &(top->idef.il[F_DISRES]);
 460     iatom   = disres->iatoms;
 461     for (i = j = k = 0; (i < disres->nr); )
 462     {
 463         type  = iatom[i];
 464         ftype = top->idef.functype[type];
 465         natom = interaction_function[ftype].nratoms+1;
 466         if (label1 != top->idef.iparams[type].disres.label)
 467         {
 468             pair[j] = k;
 469             label1  = top->idef.iparams[type].disres.label;
 470             j++;
 471         }
 472         k++;
 473         i += natom;
 474     }
 475     pair[j]  = k;
 476     *npairs  = k;
 477     if (j != nb)
 478     {
 479         gmx_incons("get_bounds for distance restraints");
 480     }
 481
 482     *index   = ind;
 483     *dr_pair = pair;
 484
 485     return nb;
 486 }
 487
 488 static void calc_violations(real rt[], real rav3[], int nb, int index[],
 489                             real bounds[], real *viol, double *st, double *sa)
 490 {
 491     const   real sixth = 1.0/6.0;
 492     int          i, j;
 493     double       rsum, rav, sumaver, sumt;
 494
 495     sumaver = 0;
 496     sumt    = 0;
 497     for (i = 0; (i < nb); i++)
 498     {
 499         rsum = 0.0;
 500         rav  = 0.0;
 501         for (j = index[i]; (j < index[i+1]); j++)
 502         {
 503             if (viol)
 504             {
 505                 viol[j] += mypow(rt[j], -3.0);
 506             }
 507             rav     += sqr(rav3[j]);
 508             rsum    += mypow(rt[j], -6);
 509         }
 510         rsum    = max(0.0, mypow(rsum, -sixth)-bounds[i]);
 511         rav     = max(0.0, mypow(rav, -sixth)-bounds[i]);
 512
 513         sumt    += rsum;
 514         sumaver += rav;
 515     }
 516     *st = sumt;
 517     *sa = sumaver;
 518 }
 519
 520 static void analyse_disre(const char *voutfn,    int nframes,
 521                           real violaver[], real bounds[], int index[],
 522                           int pair[],      int nbounds,
 523                           const output_env_t oenv)
 524 {
 525     FILE   *vout;
 526     double  sum, sumt, sumaver;
 527     int     i, j;
 528
 529     /* Subtract bounds from distances, to calculate violations */
 530     calc_violations(violaver, violaver,
 531                     nbounds, pair, bounds, NULL, &sumt, &sumaver);
 532
 533 #ifdef DEBUG
 534     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 535             sumaver);
 536     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
 537             sumt);
 538 #endif
 539     vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
 540                     oenv);
 541     sum  = 0.0;
 542     sumt = 0.0;
 543     for (i = 0; (i < nbounds); i++)
 544     {
 545         /* Do ensemble averaging */
 546         sumaver = 0;
 547         for (j = pair[i]; (j < pair[i+1]); j++)
 548         {
 549             sumaver += sqr(violaver[j]/nframes);
 550         }
 551         sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
 552
 553         sumt   += sumaver;
 554         sum     = max(sum, sumaver);
 555         fprintf(vout, "%10d  %10.5e\n", index[i], sumaver);
 556     }
 557 #ifdef DEBUG
 558     for (j = 0; (j < dr.ndr); j++)
 559     {
 560         fprintf(vout, "%10d  %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
 561     }
 562 #endif
 563     ffclose(vout);
 564
 565     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 566             sumt);
 567     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
 568
 569     do_view(oenv, voutfn, "-graphtype bar");
 570 }
 571
 572 static void einstein_visco(const char *fn, const char *fni, int nsets,
 573                            int 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     ffclose(fp0);
 626     ffclose(fp1);
 627 }
 628
 629 typedef struct {
 630     gmx_large_int_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_large_int_t nst;
 640     gmx_large_int_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_large_int_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->npoints > 0)
 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 set[], 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, gmx_bool bTempFluct,
1129                          gmx_bool bVisco, const char *visfn, int nmol,
1130                          gmx_large_int_t start_step, double start_t,
1131                          gmx_large_int_t step, double t,
1132                          double time[], 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_large_int_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->npoints == 0)
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, 1);
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, 1);
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             ffclose(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         ffclose(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_large_int_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_large_int_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_large_int) ||
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
1855         "[TT]g_energy[tt] extracts energy components or distance restraint",
1856         "data from an energy file. The user is prompted to interactively",
1857         "select the desired energy terms.[PAR]",
1858
1859         "Average, RMSD, and drift are calculated with full precision from the",
1860         "simulation (see printed manual). Drift is calculated by performing",
1861         "a least-squares fit of the data to a straight line. The reported total drift",
1862         "is the difference of the fit at the first and last point.",
1863         "An error estimate of the average is given based on a block averages",
1864         "over 5 blocks using the full-precision averages. The error estimate",
1865         "can be performed over multiple block lengths with the options",
1866         "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1867         "[BB]Note[bb] that in most cases the energy files contains averages over all",
1868         "MD steps, or over many more points than the number of frames in",
1869         "energy file. This makes the [TT]g_energy[tt] statistics output more accurate",
1870         "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1871         "file, the statistics mentioned above are simply over the single, per-frame",
1872         "energy values.[PAR]",
1873
1874         "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1875
1876         "Some fluctuation-dependent properties can be calculated provided",
1877         "the correct energy terms are selected, and that the command line option",
1878         "[TT]-fluct_props[tt] is given. The following properties",
1879         "will be computed:[BR]",
1880         "Property                        Energy terms needed[BR]",
1881         "---------------------------------------------------[BR]",
1882         "Heat capacity C[SUB]p[sub] (NPT sims):    Enthalpy, Temp     [BR]",
1883         "Heat capacity C[SUB]v[sub] (NVT sims):    Etot, Temp         [BR]",
1884         "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1885         "Isothermal compressibility:     Vol, Temp          [BR]",
1886         "Adiabatic bulk modulus:         Vol, Temp          [BR]",
1887         "---------------------------------------------------[BR]",
1888         "You always need to set the number of molecules [TT]-nmol[tt].",
1889         "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1890         "for quantum effects. Use the [TT]g_dos[tt] program if you need that (and you do).[PAR]"
1891         "When the [TT]-viol[tt] option is set, the time averaged",
1892         "violations are plotted and the running time-averaged and",
1893         "instantaneous sum of violations are recalculated. Additionally",
1894         "running time-averaged and instantaneous distances between",
1895         "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1896
1897         "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1898         "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1899         "The first two options plot the orientation, the last three the",
1900         "deviations of the orientations from the experimental values.",
1901         "The options that end on an 'a' plot the average over time",
1902         "as a function of restraint. The options that end on a 't'",
1903         "prompt the user for restraint label numbers and plot the data",
1904         "as a function of time. Option [TT]-odr[tt] plots the RMS",
1905         "deviation as a function of restraint.",
1906         "When the run used time or ensemble averaged orientation restraints,",
1907         "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1908         "not ensemble-averaged orientations and deviations instead of",
1909         "the time and ensemble averages.[PAR]",
1910
1911         "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1912         "tensor for each orientation restraint experiment. With option",
1913         "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1914
1915         "Option [TT]-odh[tt] extracts and plots the free energy data",
1916         "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1917         "from the [TT]ener.edr[tt] file.[PAR]",
1918
1919         "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1920         "difference with an ideal gas state: [BR]",
1921         "  [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
1922         "  [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
1923         "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1924         "the average is over the ensemble (or time in a trajectory).",
1925         "Note that this is in principle",
1926         "only correct when averaging over the whole (Boltzmann) ensemble",
1927         "and using the potential energy. This also allows for an entropy",
1928         "estimate using:[BR]",
1929         "  [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T[BR]",
1930         "  [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",
1931         "[PAR]",
1932
1933         "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1934         "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1935         "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1936         "files, and the average is over the ensemble A. The running average",
1937         "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1938         "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1939
1940     };
1941     static gmx_bool    bSum    = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1942     static gmx_bool    bDp     = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1943     static int         skip    = 0, nmol = 1, nbmin = 5, nbmax = 5;
1944     static real        reftemp = 300.0, ezero = 0;
1945     t_pargs            pa[]    = {
1946         { "-fee",   FALSE, etBOOL,  {&bFee},
1947           "Do a free energy estimate" },
1948         { "-fetemp", FALSE, etREAL, {&reftemp},
1949           "Reference temperature for free energy calculation" },
1950         { "-zero", FALSE, etREAL, {&ezero},
1951           "Subtract a zero-point energy" },
1952         { "-sum",  FALSE, etBOOL, {&bSum},
1953           "Sum the energy terms selected rather than display them all" },
1954         { "-dp",   FALSE, etBOOL, {&bDp},
1955           "Print energies in high precision" },
1956         { "-nbmin", FALSE, etINT, {&nbmin},
1957           "Minimum number of blocks for error estimate" },
1958         { "-nbmax", FALSE, etINT, {&nbmax},
1959           "Maximum number of blocks for error estimate" },
1960         { "-mutot", FALSE, etBOOL, {&bMutot},
1961           "Compute the total dipole moment from the components" },
1962         { "-skip", FALSE, etINT,  {&skip},
1963           "Skip number of frames between data points" },
1964         { "-aver", FALSE, etBOOL, {&bPrAll},
1965           "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1966         { "-nmol", FALSE, etINT,  {&nmol},
1967           "Number of molecules in your sample: the energies are divided by this number" },
1968         { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1969           "Compute properties based on energy fluctuations, like heat capacity" },
1970         { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1971           "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1972         { "-fluc", FALSE, etBOOL, {&bFluct},
1973           "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1974         { "-orinst", FALSE, etBOOL, {&bOrinst},
1975           "Analyse instantaneous orientation data" },
1976         { "-ovec", FALSE, etBOOL, {&bOvec},
1977           "Also plot the eigenvectors with [TT]-oten[tt]" }
1978     };
1979     const char       * drleg[] = {
1980         "Running average",
1981         "Instantaneous"
1982     };
1983     static const char *setnm[] = {
1984         "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1985         "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1986         "Volume",  "Pressure"
1987     };
1988
1989     FILE              *out     = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1990     FILE              *fp_dhdl = NULL;
1991     FILE             **drout;
1992     ener_file_t        fp;
1993     int                timecheck = 0;
1994     gmx_mtop_t         mtop;
1995     gmx_localtop_t    *top = NULL;
1996     t_inputrec         ir;
1997     t_energy         **ee;
1998     enerdata_t         edat;
1999     gmx_enxnm_t       *enm = NULL;
2000     t_enxframe        *frame, *fr = NULL;
2001     int                cur = 0;
2002 #define NEXT (1-cur)
2003     int                nre, teller, teller_disre, nfr;
2004     gmx_large_int_t    start_step;
2005     int                nor = 0, nex = 0, norfr = 0, enx_i = 0;
2006     real               start_t;
2007     real              *bounds  = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
2008     int               *index   = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
2009     int                nbounds = 0, npairs;
2010     gmx_bool           bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
2011     gmx_bool           bFoundStart, bCont, bEDR, bVisco;
2012     double             sum, sumaver, sumt, ener, dbl;
2013     double            *time = NULL;
2014     real               Vaver;
2015     int               *set     = NULL, i, j, k, nset, sss;
2016     gmx_bool          *bIsEner = NULL;
2017     char             **pairleg, **odtleg, **otenleg;
2018     char             **leg = NULL;
2019     char             **nms;
2020     char              *anm_j, *anm_k, *resnm_j, *resnm_k;
2021     int                resnr_j, resnr_k;
2022     const char        *orinst_sub = "@ subtitle \"instantaneous\"\n";
2023     char               buf[256];
2024     output_env_t       oenv;
2025     t_enxblock        *blk       = NULL;
2026     t_enxblock        *blk_disre = NULL;
2027     int                ndisre    = 0;
2028     int                dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2029
2030     t_filenm           fnm[] = {
2031         { efEDR, "-f",    NULL,      ffREAD  },
2032         { efEDR, "-f2",   NULL,      ffOPTRD },
2033         { efTPX, "-s",    NULL,      ffOPTRD },
2034         { efXVG, "-o",    "energy",  ffWRITE },
2035         { efXVG, "-viol", "violaver", ffOPTWR },
2036         { efXVG, "-pairs", "pairs",   ffOPTWR },
2037         { efXVG, "-ora",  "orienta", ffOPTWR },
2038         { efXVG, "-ort",  "orientt", ffOPTWR },
2039         { efXVG, "-oda",  "orideva", ffOPTWR },
2040         { efXVG, "-odr",  "oridevr", ffOPTWR },
2041         { efXVG, "-odt",  "oridevt", ffOPTWR },
2042         { efXVG, "-oten", "oriten",  ffOPTWR },
2043         { efXVG, "-corr", "enecorr", ffOPTWR },
2044         { efXVG, "-vis",  "visco",   ffOPTWR },
2045         { efXVG, "-ravg", "runavgdf", ffOPTWR },
2046         { efXVG, "-odh",  "dhdl", ffOPTWR }
2047     };
2048 #define NFILE asize(fnm)
2049     int                npargs;
2050     t_pargs           *ppa;
2051
2052     CopyRight(stderr, argv[0]);
2053     npargs = asize(pa);
2054     ppa    = add_acf_pargs(&npargs, pa);
2055     parse_common_args(&argc, argv,
2056                       PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
2057                       NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv);
2058
2059     bDRAll = opt2bSet("-pairs", NFILE, fnm);
2060     bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2061     bORA   = opt2bSet("-ora", NFILE, fnm);
2062     bORT   = opt2bSet("-ort", NFILE, fnm);
2063     bODA   = opt2bSet("-oda", NFILE, fnm);
2064     bODR   = opt2bSet("-odr", NFILE, fnm);
2065     bODT   = opt2bSet("-odt", NFILE, fnm);
2066     bORIRE = bORA || bORT || bODA || bODR || bODT;
2067     bOTEN  = opt2bSet("-oten", NFILE, fnm);
2068     bDHDL  = opt2bSet("-odh", NFILE, fnm);
2069
2070     nset = 0;
2071
2072     snew(frame, 2);
2073     fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2074     do_enxnms(fp, &nre, &enm);
2075
2076     Vaver = -1;
2077
2078     bVisco = opt2bSet("-vis", NFILE, fnm);
2079
2080     if ((!bDisRe) && (!bDHDL))
2081     {
2082         if (bVisco)
2083         {
2084             nset = asize(setnm);
2085             snew(set, nset);
2086             /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2087             for (j = 0; j < nset; j++)
2088             {
2089                 for (i = 0; i < nre; i++)
2090                 {
2091                     if (strstr(enm[i].name, setnm[j]))
2092                     {
2093                         set[j] = i;
2094                         break;
2095                     }
2096                 }
2097                 if (i == nre)
2098                 {
2099                     if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2100                     {
2101                         printf("Enter the box volume (" unit_volume "): ");
2102                         if (1 != scanf("%lf", &dbl))
2103                         {
2104                             gmx_fatal(FARGS, "Error reading user input");
2105                         }
2106                         Vaver = dbl;
2107                     }
2108                     else
2109                     {
2110                         gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2111                                   setnm[j]);
2112                     }
2113                 }
2114             }
2115         }
2116         else
2117         {
2118             set = select_by_name(nre, enm, &nset);
2119         }
2120         /* Print all the different units once */
2121         sprintf(buf, "(%s)", enm[set[0]].unit);
2122         for (i = 1; i < nset; i++)
2123         {
2124             for (j = 0; j < i; j++)
2125             {
2126                 if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2127                 {
2128                     break;
2129                 }
2130             }
2131             if (j == i)
2132             {
2133                 strcat(buf, ", (");
2134                 strcat(buf, enm[set[i]].unit);
2135                 strcat(buf, ")");
2136             }
2137         }
2138         out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf,
2139                        oenv);
2140
2141         snew(leg, nset+1);
2142         for (i = 0; (i < nset); i++)
2143         {
2144             leg[i] = enm[set[i]].name;
2145         }
2146         if (bSum)
2147         {
2148             leg[nset] = strdup("Sum");
2149             xvgr_legend(out, nset+1, (const char**)leg, oenv);
2150         }
2151         else
2152         {
2153             xvgr_legend(out, nset, (const char**)leg, oenv);
2154         }
2155
2156         snew(bIsEner, nset);
2157         for (i = 0; (i < nset); i++)
2158         {
2159             bIsEner[i] = FALSE;
2160             for (j = 0; (j <= F_ETOT); j++)
2161             {
2162                 bIsEner[i] = bIsEner[i] ||
2163                     (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2164             }
2165         }
2166
2167         if (bPrAll && nset > 1)
2168         {
2169             gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2170         }
2171
2172         time = NULL;
2173
2174         if (bORIRE || bOTEN)
2175         {
2176             get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2177         }
2178
2179         if (bORIRE)
2180         {
2181             if (bOrinst)
2182             {
2183                 enx_i = enxORI;
2184             }
2185             else
2186             {
2187                 enx_i = enxOR;
2188             }
2189
2190             if (bORA || bODA)
2191             {
2192                 snew(orient, nor);
2193             }
2194             if (bODR)
2195             {
2196                 snew(odrms, nor);
2197             }
2198             if (bORT || bODT)
2199             {
2200                 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2201                 fprintf(stderr, "End your selection with 0\n");
2202                 j     = -1;
2203                 orsel = NULL;
2204                 do
2205                 {
2206                     j++;
2207                     srenew(orsel, j+1);
2208                     if (1 != scanf("%d", &(orsel[j])))
2209                     {
2210                         gmx_fatal(FARGS, "Error reading user input");
2211                     }
2212                 }
2213                 while (orsel[j] > 0);
2214                 if (orsel[0] == -1)
2215                 {
2216                     fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2217                     norsel = nor;
2218                     srenew(orsel, nor);
2219                     for (i = 0; i < nor; i++)
2220                     {
2221                         orsel[i] = i;
2222                     }
2223                 }
2224                 else
2225                 {
2226                     /* Build the selection */
2227                     norsel = 0;
2228                     for (i = 0; i < j; i++)
2229                     {
2230                         for (k = 0; k < nor; k++)
2231                         {
2232                             if (or_label[k] == orsel[i])
2233                             {
2234                                 orsel[norsel] = k;
2235                                 norsel++;
2236                                 break;
2237                             }
2238                         }
2239                         if (k == nor)
2240                         {
2241                             fprintf(stderr, "Orientation restraint label %d not found\n",
2242                                     orsel[i]);
2243                         }
2244                     }
2245                 }
2246                 snew(odtleg, norsel);
2247                 for (i = 0; i < norsel; i++)
2248                 {
2249                     snew(odtleg[i], 256);
2250                     sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2251                 }
2252                 if (bORT)
2253                 {
2254                     fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2255                                     "Time (ps)", "", oenv);
2256                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2257                     {
2258                         fprintf(fort, "%s", orinst_sub);
2259                     }
2260                     xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2261                 }
2262                 if (bODT)
2263                 {
2264                     fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2265                                     "Orientation restraint deviation",
2266                                     "Time (ps)", "", oenv);
2267                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2268                     {
2269                         fprintf(fodt, "%s", orinst_sub);
2270                     }
2271                     xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2272                 }
2273             }
2274         }
2275         if (bOTEN)
2276         {
2277             foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2278                              "Order tensor", "Time (ps)", "", oenv);
2279             snew(otenleg, bOvec ? nex*12 : nex*3);
2280             for (i = 0; i < nex; i++)
2281             {
2282                 for (j = 0; j < 3; j++)
2283                 {
2284                     sprintf(buf, "eig%d", j+1);
2285                     otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
2286                 }
2287                 if (bOvec)
2288                 {
2289                     for (j = 0; j < 9; j++)
2290                     {
2291                         sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2292                         otenleg[12*i+3+j] = strdup(buf);
2293                     }
2294                 }
2295             }
2296             xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2297         }
2298     }
2299     else if (bDisRe)
2300     {
2301         nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
2302                              &mtop, &top, &ir);
2303         snew(violaver, npairs);
2304         out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2305                        "Time (ps)", "nm", oenv);
2306         xvgr_legend(out, 2, drleg, oenv);
2307         if (bDRAll)
2308         {
2309             fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2310                                 "Time (ps)", "Distance (nm)", oenv);
2311             if (output_env_get_print_xvgr_codes(oenv))
2312             {
2313                 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2314                         ir.dr_tau);
2315             }
2316         }
2317     }
2318     else if (bDHDL)
2319     {
2320         get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
2321     }
2322
2323     /* Initiate energies and set them to zero */
2324     edat.nsteps  = 0;
2325     edat.npoints = 0;
2326     edat.nframes = 0;
2327     edat.step    = NULL;
2328     edat.steps   = NULL;
2329     edat.points  = NULL;
2330     snew(edat.s, nset);
2331
2332     /* Initiate counters */
2333     teller       = 0;
2334     teller_disre = 0;
2335     bFoundStart  = FALSE;
2336     start_step   = 0;
2337     start_t      = 0;
2338     do
2339     {
2340         /* This loop searches for the first frame (when -b option is given),
2341          * or when this has been found it reads just one energy frame
2342          */
2343         do
2344         {
2345             bCont = do_enx(fp, &(frame[NEXT]));
2346             if (bCont)
2347             {
2348                 timecheck = check_times(frame[NEXT].t);
2349             }
2350         }
2351         while (bCont && (timecheck < 0));
2352
2353         if ((timecheck == 0) && bCont)
2354         {
2355             /* We read a valid frame, so we can use it */
2356             fr = &(frame[NEXT]);
2357
2358             if (fr->nre > 0)
2359             {
2360                 /* The frame contains energies, so update cur */
2361                 cur  = NEXT;
2362
2363                 if (edat.nframes % 1000 == 0)
2364                 {
2365                     srenew(edat.step, edat.nframes+1000);
2366                     memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2367                     srenew(edat.steps, edat.nframes+1000);
2368                     memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2369                     srenew(edat.points, edat.nframes+1000);
2370                     memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2371
2372                     for (i = 0; i < nset; i++)
2373                     {
2374                         srenew(edat.s[i].ener, edat.nframes+1000);
2375                         memset(&(edat.s[i].ener[edat.nframes]), 0,
2376                                1000*sizeof(edat.s[i].ener[0]));
2377                         srenew(edat.s[i].es, edat.nframes+1000);
2378                         memset(&(edat.s[i].es[edat.nframes]), 0,
2379                                1000*sizeof(edat.s[i].es[0]));
2380                     }
2381                 }
2382
2383                 nfr            = edat.nframes;
2384                 edat.step[nfr] = fr->step;
2385
2386                 if (!bFoundStart)
2387                 {
2388                     bFoundStart = TRUE;
2389                     /* Initiate the previous step data */
2390                     start_step = fr->step;
2391                     start_t    = fr->t;
2392                     /* Initiate the energy sums */
2393                     edat.steps[nfr]  = 1;
2394                     edat.points[nfr] = 1;
2395                     for (i = 0; i < nset; i++)
2396                     {
2397                         sss                    = set[i];
2398                         edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2399                         edat.s[i].es[nfr].sum2 = 0;
2400                     }
2401                     edat.nsteps  = 1;
2402                     edat.npoints = 1;
2403                 }
2404                 else
2405                 {
2406                     edat.steps[nfr] = fr->nsteps;
2407                     {
2408                         if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2409                         {
2410                             if (fr->nsum <= 1)
2411                             {
2412                                 edat.points[nfr] = 1;
2413                                 for (i = 0; i < nset; i++)
2414                                 {
2415                                     sss                    = set[i];
2416                                     edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2417                                     edat.s[i].es[nfr].sum2 = 0;
2418                                 }
2419                                 edat.npoints += 1;
2420                             }
2421                             else
2422                             {
2423                                 edat.points[nfr] = fr->nsum;
2424                                 for (i = 0; i < nset; i++)
2425                                 {
2426                                     sss                    = set[i];
2427                                     edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2428                                     edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2429                                 }
2430                                 edat.npoints += fr->nsum;
2431                             }
2432                         }
2433                         else
2434                         {
2435                             /* The interval does not match fr->nsteps:
2436                              * can not do exact averages.
2437                              */
2438                             edat.npoints = 0;
2439                         }
2440                         edat.nsteps = fr->step - start_step + 1;
2441                     }
2442                 }
2443                 for (i = 0; i < nset; i++)
2444                 {
2445                     edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2446                 }
2447             }
2448             /*
2449              * Define distance restraint legends. Can only be done after
2450              * the first frame has been read... (Then we know how many there are)
2451              */
2452             blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2453             if (bDisRe && bDRAll && !leg && blk_disre)
2454             {
2455                 t_iatom   *fa;
2456                 t_iparams *ip;
2457
2458                 fa = top->idef.il[F_DISRES].iatoms;
2459                 ip = top->idef.iparams;
2460                 if (blk_disre->nsub != 2 ||
2461                     (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2462                 {
2463                     gmx_incons("Number of disre sub-blocks not equal to 2");
2464                 }
2465
2466                 ndisre = blk_disre->sub[0].nr;
2467                 if (ndisre != top->idef.il[F_DISRES].nr/3)
2468                 {
2469                     gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2470                               ndisre, top->idef.il[F_DISRES].nr/3);
2471                 }
2472                 snew(pairleg, ndisre);
2473                 for (i = 0; i < ndisre; i++)
2474                 {
2475                     snew(pairleg[i], 30);
2476                     j = fa[3*i+1];
2477                     k = fa[3*i+2];
2478                     gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2479                     gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2480                     sprintf(pairleg[i], "%d %s %d %s (%d)",
2481                             resnr_j, anm_j, resnr_k, anm_k,
2482                             ip[fa[3*i]].disres.label);
2483                 }
2484                 set = select_it(ndisre, pairleg, &nset);
2485                 snew(leg, 2*nset);
2486                 for (i = 0; (i < nset); i++)
2487                 {
2488                     snew(leg[2*i], 32);
2489                     sprintf(leg[2*i],  "a %s", pairleg[set[i]]);
2490                     snew(leg[2*i+1], 32);
2491                     sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2492                 }
2493                 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2494             }
2495
2496             /*
2497              * Store energies for analysis afterwards...
2498              */
2499             if (!bDisRe && !bDHDL && (fr->nre > 0))
2500             {
2501                 if (edat.nframes % 1000 == 0)
2502                 {
2503                     srenew(time, edat.nframes+1000);
2504                 }
2505                 time[edat.nframes] = fr->t;
2506                 edat.nframes++;
2507             }
2508             /*
2509              * Printing time, only when we do not want to skip frames
2510              */
2511             if (!skip || teller % skip == 0)
2512             {
2513                 if (bDisRe)
2514                 {
2515                     /*******************************************
2516                      * D I S T A N C E   R E S T R A I N T S
2517                      *******************************************/
2518                     if (ndisre > 0)
2519                     {
2520  #ifndef GMX_DOUBLE
2521                         float  *disre_rt     =     blk_disre->sub[0].fval;
2522                         float  *disre_rm3tav = blk_disre->sub[1].fval;
2523  #else
2524                         double *disre_rt     =     blk_disre->sub[0].dval;
2525                         double *disre_rm3tav = blk_disre->sub[1].dval;
2526  #endif
2527
2528                         print_time(out, fr->t);
2529                         if (violaver == NULL)
2530                         {
2531                             snew(violaver, ndisre);
2532                         }
2533
2534                         /* Subtract bounds from distances, to calculate violations */
2535                         calc_violations(disre_rt, disre_rm3tav,
2536                                         nbounds, pair, bounds, violaver, &sumt, &sumaver);
2537
2538                         fprintf(out, "  %8.4f  %8.4f\n", sumaver, sumt);
2539                         if (bDRAll)
2540                         {
2541                             print_time(fp_pairs, fr->t);
2542                             for (i = 0; (i < nset); i++)
2543                             {
2544                                 sss = set[i];
2545                                 fprintf(fp_pairs, "  %8.4f", mypow(disre_rm3tav[sss], minthird));
2546                                 fprintf(fp_pairs, "  %8.4f", disre_rt[sss]);
2547                             }
2548                             fprintf(fp_pairs, "\n");
2549                         }
2550                         teller_disre++;
2551                     }
2552                 }
2553                 else if (bDHDL)
2554                 {
2555                     do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2556                 }
2557
2558                 /*******************************************
2559                  * E N E R G I E S
2560                  *******************************************/
2561                 else
2562                 {
2563                     if (fr->nre > 0)
2564                     {
2565                         if (bPrAll)
2566                         {
2567                             /* We skip frames with single points (usually only the first frame),
2568                              * since they would result in an average plot with outliers.
2569                              */
2570                             if (fr->nsum > 1)
2571                             {
2572                                 print_time(out, fr->t);
2573                                 print1(out, bDp, fr->ener[set[0]].e);
2574                                 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2575                                 print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
2576                                 fprintf(out, "\n");
2577                             }
2578                         }
2579                         else
2580                         {
2581                             print_time(out, fr->t);
2582                             if (bSum)
2583                             {
2584                                 sum = 0;
2585                                 for (i = 0; i < nset; i++)
2586                                 {
2587                                     sum += fr->ener[set[i]].e;
2588                                 }
2589                                 print1(out, bDp, sum/nmol-ezero);
2590                             }
2591                             else
2592                             {
2593                                 for (i = 0; (i < nset); i++)
2594                                 {
2595                                     if (bIsEner[i])
2596                                     {
2597                                         print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2598                                     }
2599                                     else
2600                                     {
2601                                         print1(out, bDp, fr->ener[set[i]].e);
2602                                     }
2603                                 }
2604                             }
2605                             fprintf(out, "\n");
2606                         }
2607                     }
2608                     blk = find_block_id_enxframe(fr, enx_i, NULL);
2609                     if (bORIRE && blk)
2610                     {
2611 #ifndef GMX_DOUBLE
2612                         xdr_datatype dt = xdr_datatype_float;
2613 #else
2614                         xdr_datatype dt = xdr_datatype_double;
2615 #endif
2616                         real        *vals;
2617
2618                         if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2619                         {
2620                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2621                         }
2622 #ifndef GMX_DOUBLE
2623                         vals = blk->sub[0].fval;
2624 #else
2625                         vals = blk->sub[0].dval;
2626 #endif
2627
2628                         if (blk->sub[0].nr != (size_t)nor)
2629                         {
2630                             gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2631                         }
2632                         if (bORA || bODA)
2633                         {
2634                             for (i = 0; i < nor; i++)
2635                             {
2636                                 orient[i] += vals[i];
2637                             }
2638                         }
2639                         if (bODR)
2640                         {
2641                             for (i = 0; i < nor; i++)
2642                             {
2643                                 odrms[i] += sqr(vals[i]-oobs[i]);
2644                             }
2645                         }
2646                         if (bORT)
2647                         {
2648                             fprintf(fort, "  %10f", fr->t);
2649                             for (i = 0; i < norsel; i++)
2650                             {
2651                                 fprintf(fort, " %g", vals[orsel[i]]);
2652                             }
2653                             fprintf(fort, "\n");
2654                         }
2655                         if (bODT)
2656                         {
2657                             fprintf(fodt, "  %10f", fr->t);
2658                             for (i = 0; i < norsel; i++)
2659                             {
2660                                 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2661                             }
2662                             fprintf(fodt, "\n");
2663                         }
2664                         norfr++;
2665                     }
2666                     blk = find_block_id_enxframe(fr, enxORT, NULL);
2667                     if (bOTEN && blk)
2668                     {
2669 #ifndef GMX_DOUBLE
2670                         xdr_datatype dt = xdr_datatype_float;
2671 #else
2672                         xdr_datatype dt = xdr_datatype_double;
2673 #endif
2674                         real        *vals;
2675
2676                         if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2677                         {
2678                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2679                         }
2680 #ifndef GMX_DOUBLE
2681                         vals = blk->sub[0].fval;
2682 #else
2683                         vals = blk->sub[0].dval;
2684 #endif
2685
2686                         if (blk->sub[0].nr != (size_t)(nex*12))
2687                         {
2688                             gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2689                                       blk->sub[0].nr/12, nex);
2690                         }
2691                         fprintf(foten, "  %10f", fr->t);
2692                         for (i = 0; i < nex; i++)
2693                         {
2694                             for (j = 0; j < (bOvec ? 12 : 3); j++)
2695                             {
2696                                 fprintf(foten, " %g", vals[i*12+j]);
2697                             }
2698                         }
2699                         fprintf(foten, "\n");
2700                     }
2701                 }
2702             }
2703             teller++;
2704         }
2705     }
2706     while (bCont && (timecheck == 0));
2707
2708     fprintf(stderr, "\n");
2709     close_enx(fp);
2710     if (out)
2711     {
2712         ffclose(out);
2713     }
2714
2715     if (bDRAll)
2716     {
2717         ffclose(fp_pairs);
2718     }
2719
2720     if (bORT)
2721     {
2722         ffclose(fort);
2723     }
2724     if (bODT)
2725     {
2726         ffclose(fodt);
2727     }
2728     if (bORA)
2729     {
2730         out = xvgropen(opt2fn("-ora", NFILE, fnm),
2731                        "Average calculated orientations",
2732                        "Restraint label", "", oenv);
2733         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2734         {
2735             fprintf(out, "%s", orinst_sub);
2736         }
2737         for (i = 0; i < nor; i++)
2738         {
2739             fprintf(out, "%5d  %g\n", or_label[i], orient[i]/norfr);
2740         }
2741         ffclose(out);
2742     }
2743     if (bODA)
2744     {
2745         out = xvgropen(opt2fn("-oda", NFILE, fnm),
2746                        "Average restraint deviation",
2747                        "Restraint label", "", oenv);
2748         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2749         {
2750             fprintf(out, "%s", orinst_sub);
2751         }
2752         for (i = 0; i < nor; i++)
2753         {
2754             fprintf(out, "%5d  %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2755         }
2756         ffclose(out);
2757     }
2758     if (bODR)
2759     {
2760         out = xvgropen(opt2fn("-odr", NFILE, fnm),
2761                        "RMS orientation restraint deviations",
2762                        "Restraint label", "", oenv);
2763         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2764         {
2765             fprintf(out, "%s", orinst_sub);
2766         }
2767         for (i = 0; i < nor; i++)
2768         {
2769             fprintf(out, "%5d  %g\n", or_label[i], sqrt(odrms[i]/norfr));
2770         }
2771         ffclose(out);
2772     }
2773     if (bOTEN)
2774     {
2775         ffclose(foten);
2776     }
2777
2778     if (bDisRe)
2779     {
2780         analyse_disre(opt2fn("-viol", NFILE, fnm),
2781                       teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2782     }
2783     else if (bDHDL)
2784     {
2785         if (fp_dhdl)
2786         {
2787             ffclose(fp_dhdl);
2788             printf("\n\nWrote %d lambda values with %d samples as ",
2789                    dh_lambdas, dh_samples);
2790             if (dh_hists > 0)
2791             {
2792                 printf("%d dH histograms ", dh_hists);
2793             }
2794             if (dh_blocks > 0)
2795             {
2796                 printf("%d dH data blocks ", dh_blocks);
2797             }
2798             printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2799
2800         }
2801         else
2802         {
2803             gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2804         }
2805
2806     }
2807     else
2808     {
2809         double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2810         analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2811                      bFee, bSum, bFluct, opt2parg_bSet("-nmol", npargs, ppa),
2812                      bVisco, opt2fn("-vis", NFILE, fnm),
2813                      nmol,
2814                      start_step, start_t, frame[cur].step, frame[cur].t,
2815                      time, reftemp, &edat,
2816                      nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2817                      oenv);
2818         if (bFluctProps)
2819         {
2820             calc_fluctuation_props(stdout, bDriftCorr, dt, nset, set, nmol, leg, &edat,
2821                                    nbmin, nbmax);
2822         }
2823     }
2824     if (opt2bSet("-f2", NFILE, fnm))
2825     {
2826         fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2827             reftemp, nset, set, leg, &edat, time, oenv);
2828     }
2829
2830     {
2831         const char *nxy = "-nxy";
2832
2833         do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2834         do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2835         do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2836         do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2837         do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2838         do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2839         do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2840         do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2841         do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);
2842     }
2843     thanx(stderr);
2844
2845     return 0;
2846 }