src/gromacs/gmxana/gmx_energy.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <math.h>
  42 #include <stdlib.h>
  43 #include <string.h>
  44
  45 #include "gromacs/legacyheaders/typedefs.h"
  46 #include "gromacs/utility/fatalerror.h"
  47 #include "gromacs/math/vec.h"
  48 #include "gromacs/utility/cstringutil.h"
  49 #include "gromacs/utility/smalloc.h"
  50 #include "gromacs/fileio/enxio.h"
  51 #include "gromacs/commandline/pargs.h"
  52 #include "gromacs/legacyheaders/names.h"
  53 #include "gromacs/legacyheaders/copyrite.h"
  54 #include "gromacs/legacyheaders/macros.h"
  55 #include "gromacs/fileio/xvgr.h"
  56 #include "gstat.h"
  57 #include "gromacs/math/units.h"
  58 #include "gromacs/fileio/tpxio.h"
  59 #include "gromacs/fileio/trxio.h"
  60 #include "gromacs/legacyheaders/viewit.h"
  61 #include "gromacs/topology/mtop_util.h"
  62 #include "gmx_ana.h"
  63 #include "gromacs/legacyheaders/mdebin.h"
  64
  65 static real       minthird = -1.0/3.0, minsixth = -1.0/6.0;
  66
  67 typedef struct {
  68     real sum;
  69     real sum2;
  70 } exactsum_t;
  71
  72 typedef struct {
  73     real       *ener;
  74     exactsum_t *es;
  75     gmx_bool    bExactStat;
  76     double      av;
  77     double      rmsd;
  78     double      ee;
  79     double      slope;
  80 } enerdat_t;
  81
  82 typedef struct {
  83     gmx_int64_t      nsteps;
  84     gmx_int64_t      npoints;
  85     int              nframes;
  86     int             *step;
  87     int             *steps;
  88     int             *points;
  89     enerdat_t       *s;
  90 } enerdata_t;
  91
  92 static double mypow(double x, double y)
  93 {
  94     if (x > 0)
  95     {
  96         return pow(x, y);
  97     }
  98     else
  99     {
 100         return 0.0;
 101     }
 102 }
 103
 104 static int *select_it(int nre, char *nm[], int *nset)
 105 {
 106     gmx_bool *bE;
 107     int       n, k, j, i;
 108     int      *set;
 109     gmx_bool  bVerbose = TRUE;
 110
 111     if ((getenv("GMX_ENER_VERBOSE")) != NULL)
 112     {
 113         bVerbose = FALSE;
 114     }
 115
 116     fprintf(stderr, "Select the terms you want from the following list\n");
 117     fprintf(stderr, "End your selection with 0\n");
 118
 119     if (bVerbose)
 120     {
 121         for (k = 0; (k < nre); )
 122         {
 123             for (j = 0; (j < 4) && (k < nre); j++, k++)
 124             {
 125                 fprintf(stderr, " %3d=%14s", k+1, nm[k]);
 126             }
 127             fprintf(stderr, "\n");
 128         }
 129     }
 130
 131     snew(bE, nre);
 132     do
 133     {
 134         if (1 != scanf("%d", &n))
 135         {
 136             gmx_fatal(FARGS, "Error reading user input");
 137         }
 138         if ((n > 0) && (n <= nre))
 139         {
 140             bE[n-1] = TRUE;
 141         }
 142     }
 143     while (n != 0);
 144
 145     snew(set, nre);
 146     for (i = (*nset) = 0; (i < nre); i++)
 147     {
 148         if (bE[i])
 149         {
 150             set[(*nset)++] = i;
 151         }
 152     }
 153
 154     sfree(bE);
 155
 156     return set;
 157 }
 158
 159 static void chomp(char *buf)
 160 {
 161     int len = strlen(buf);
 162
 163     while ((len > 0) && (buf[len-1] == '\n'))
 164     {
 165         buf[len-1] = '\0';
 166         len--;
 167     }
 168 }
 169
 170 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
 171 {
 172     gmx_bool   *bE;
 173     int         n, k, kk, j, i, nmatch, nind, nss;
 174     int        *set;
 175     gmx_bool    bEOF, bVerbose = TRUE, bLong = FALSE;
 176     char       *ptr, buf[STRLEN];
 177     const char *fm4   = "%3d  %-14s";
 178     const char *fm2   = "%3d  %-34s";
 179     char      **newnm = NULL;
 180
 181     if ((getenv("GMX_ENER_VERBOSE")) != NULL)
 182     {
 183         bVerbose = FALSE;
 184     }
 185
 186     fprintf(stderr, "\n");
 187     fprintf(stderr, "Select the terms you want from the following list by\n");
 188     fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
 189     fprintf(stderr, "End your selection with an empty line or a zero.\n");
 190     fprintf(stderr, "-------------------------------------------------------------------\n");
 191
 192     snew(newnm, nre);
 193     j = 0;
 194     for (k = 0; k < nre; k++)
 195     {
 196         newnm[k] = gmx_strdup(nm[k].name);
 197         /* Insert dashes in all the names */
 198         while ((ptr = strchr(newnm[k], ' ')) != NULL)
 199         {
 200             *ptr = '-';
 201         }
 202         if (bVerbose)
 203         {
 204             if (j == 0)
 205             {
 206                 if (k > 0)
 207                 {
 208                     fprintf(stderr, "\n");
 209                 }
 210                 bLong = FALSE;
 211                 for (kk = k; kk < k+4; kk++)
 212                 {
 213                     if (kk < nre && strlen(nm[kk].name) > 14)
 214                     {
 215                         bLong = TRUE;
 216                     }
 217                 }
 218             }
 219             else
 220             {
 221                 fprintf(stderr, " ");
 222             }
 223             if (!bLong)
 224             {
 225                 fprintf(stderr, fm4, k+1, newnm[k]);
 226                 j++;
 227                 if (j == 4)
 228                 {
 229                     j = 0;
 230                 }
 231             }
 232             else
 233             {
 234                 fprintf(stderr, fm2, k+1, newnm[k]);
 235                 j++;
 236                 if (j == 2)
 237                 {
 238                     j = 0;
 239                 }
 240             }
 241         }
 242     }
 243     if (bVerbose)
 244     {
 245         fprintf(stderr, "\n\n");
 246     }
 247
 248     snew(bE, nre);
 249
 250     bEOF = FALSE;
 251     while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
 252     {
 253         /* Remove newlines */
 254         chomp(buf);
 255
 256         /* Remove spaces */
 257         trim(buf);
 258
 259         /* Empty line means end of input */
 260         bEOF = (strlen(buf) == 0);
 261         if (!bEOF)
 262         {
 263             ptr = buf;
 264             do
 265             {
 266                 if (!bEOF)
 267                 {
 268                     /* First try to read an integer */
 269                     nss   = sscanf(ptr, "%d", &nind);
 270                     if (nss == 1)
 271                     {
 272                         /* Zero means end of input */
 273                         if (nind == 0)
 274                         {
 275                             bEOF = TRUE;
 276                         }
 277                         else if ((1 <= nind) && (nind <= nre))
 278                         {
 279                             bE[nind-1] = TRUE;
 280                         }
 281                         else
 282                         {
 283                             fprintf(stderr, "number %d is out of range\n", nind);
 284                         }
 285                     }
 286                     else
 287                     {
 288                         /* Now try to read a string */
 289                         i      = strlen(ptr);
 290                         nmatch = 0;
 291                         for (nind = 0; nind < nre; nind++)
 292                         {
 293                             if (gmx_strcasecmp(newnm[nind], ptr) == 0)
 294                             {
 295                                 bE[nind] = TRUE;
 296                                 nmatch++;
 297                             }
 298                         }
 299                         if (nmatch == 0)
 300                         {
 301                             i      = strlen(ptr);
 302                             nmatch = 0;
 303                             for (nind = 0; nind < nre; nind++)
 304                             {
 305                                 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
 306                                 {
 307                                     bE[nind] = TRUE;
 308                                     nmatch++;
 309                                 }
 310                             }
 311                             if (nmatch == 0)
 312                             {
 313                                 fprintf(stderr, "String '%s' does not match anything\n", ptr);
 314                             }
 315                         }
 316                     }
 317                 }
 318                 /* Look for the first space, and remove spaces from there */
 319                 if ((ptr = strchr(ptr, ' ')) != NULL)
 320                 {
 321                     trim(ptr);
 322                 }
 323             }
 324             while (!bEOF && (ptr && (strlen(ptr) > 0)));
 325         }
 326     }
 327
 328     snew(set, nre);
 329     for (i = (*nset) = 0; (i < nre); i++)
 330     {
 331         if (bE[i])
 332         {
 333             set[(*nset)++] = i;
 334         }
 335     }
 336
 337     sfree(bE);
 338
 339     if (*nset == 0)
 340     {
 341         gmx_fatal(FARGS, "No energy terms selected");
 342     }
 343
 344     for (i = 0; (i < nre); i++)
 345     {
 346         sfree(newnm[i]);
 347     }
 348     sfree(newnm);
 349
 350     return set;
 351 }
 352
 353 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
 354 {
 355     gmx_mtop_t  mtop;
 356     int         natoms;
 357     t_iatom    *iatom;
 358     matrix      box;
 359
 360     /* all we need is the ir to be able to write the label */
 361     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, &mtop);
 362 }
 363
 364 static void get_orires_parms(const char *topnm,
 365                              int *nor, int *nex, int **label, real **obs)
 366 {
 367     gmx_mtop_t      mtop;
 368     gmx_localtop_t *top;
 369     t_inputrec      ir;
 370     t_iparams      *ip;
 371     int             natoms, i;
 372     t_iatom        *iatom;
 373     int             nb;
 374     matrix          box;
 375
 376     read_tpx(topnm, &ir, box, &natoms, NULL, NULL, NULL, &mtop);
 377     top = gmx_mtop_generate_local_top(&mtop, &ir);
 378
 379     ip       = top->idef.iparams;
 380     iatom    = top->idef.il[F_ORIRES].iatoms;
 381
 382     /* Count how many distance restraint there are... */
 383     nb = top->idef.il[F_ORIRES].nr;
 384     if (nb == 0)
 385     {
 386         gmx_fatal(FARGS, "No orientation restraints in topology!\n");
 387     }
 388
 389     *nor = nb/3;
 390     *nex = 0;
 391     snew(*label, *nor);
 392     snew(*obs, *nor);
 393     for (i = 0; i < nb; i += 3)
 394     {
 395         (*label)[i/3] = ip[iatom[i]].orires.label;
 396         (*obs)[i/3]   = ip[iatom[i]].orires.obs;
 397         if (ip[iatom[i]].orires.ex >= *nex)
 398         {
 399             *nex = ip[iatom[i]].orires.ex+1;
 400         }
 401     }
 402     fprintf(stderr, "Found %d orientation restraints with %d experiments",
 403             *nor, *nex);
 404 }
 405
 406 static int get_bounds(const char *topnm,
 407                       real **bounds, int **index, int **dr_pair, int *npairs,
 408                       gmx_mtop_t *mtop, gmx_localtop_t **ltop, t_inputrec *ir)
 409 {
 410     gmx_localtop_t *top;
 411     t_functype     *functype;
 412     t_iparams      *ip;
 413     int             natoms, i, j, k, type, ftype, natom;
 414     t_ilist        *disres;
 415     t_iatom        *iatom;
 416     real           *b;
 417     int            *ind, *pair;
 418     int             nb, label1;
 419     matrix          box;
 420
 421     read_tpx(topnm, ir, box, &natoms, NULL, NULL, NULL, mtop);
 422     snew(*ltop, 1);
 423     top   = gmx_mtop_generate_local_top(mtop, ir);
 424     *ltop = top;
 425
 426     functype = top->idef.functype;
 427     ip       = top->idef.iparams;
 428
 429     /* Count how many distance restraint there are... */
 430     nb = top->idef.il[F_DISRES].nr;
 431     if (nb == 0)
 432     {
 433         gmx_fatal(FARGS, "No distance restraints in topology!\n");
 434     }
 435
 436     /* Allocate memory */
 437     snew(b, nb);
 438     snew(ind, nb);
 439     snew(pair, nb+1);
 440
 441     /* Fill the bound array */
 442     nb = 0;
 443     for (i = 0; (i < top->idef.ntypes); i++)
 444     {
 445         ftype = functype[i];
 446         if (ftype == F_DISRES)
 447         {
 448
 449             label1  = ip[i].disres.label;
 450             b[nb]   = ip[i].disres.up1;
 451             ind[nb] = label1;
 452             nb++;
 453         }
 454     }
 455     *bounds = b;
 456
 457     /* Fill the index array */
 458     label1  = -1;
 459     disres  = &(top->idef.il[F_DISRES]);
 460     iatom   = disres->iatoms;
 461     for (i = j = k = 0; (i < disres->nr); )
 462     {
 463         type  = iatom[i];
 464         ftype = top->idef.functype[type];
 465         natom = interaction_function[ftype].nratoms+1;
 466         if (label1 != top->idef.iparams[type].disres.label)
 467         {
 468             pair[j] = k;
 469             label1  = top->idef.iparams[type].disres.label;
 470             j++;
 471         }
 472         k++;
 473         i += natom;
 474     }
 475     pair[j]  = k;
 476     *npairs  = k;
 477     if (j != nb)
 478     {
 479         gmx_incons("get_bounds for distance restraints");
 480     }
 481
 482     *index   = ind;
 483     *dr_pair = pair;
 484
 485     return nb;
 486 }
 487
 488 static void calc_violations(real rt[], real rav3[], int nb, int index[],
 489                             real bounds[], real *viol, double *st, double *sa)
 490 {
 491     const   real sixth = 1.0/6.0;
 492     int          i, j;
 493     double       rsum, rav, sumaver, sumt;
 494
 495     sumaver = 0;
 496     sumt    = 0;
 497     for (i = 0; (i < nb); i++)
 498     {
 499         rsum = 0.0;
 500         rav  = 0.0;
 501         for (j = index[i]; (j < index[i+1]); j++)
 502         {
 503             if (viol)
 504             {
 505                 viol[j] += mypow(rt[j], -3.0);
 506             }
 507             rav     += sqr(rav3[j]);
 508             rsum    += mypow(rt[j], -6);
 509         }
 510         rsum    = max(0.0, mypow(rsum, -sixth)-bounds[i]);
 511         rav     = max(0.0, mypow(rav, -sixth)-bounds[i]);
 512
 513         sumt    += rsum;
 514         sumaver += rav;
 515     }
 516     *st = sumt;
 517     *sa = sumaver;
 518 }
 519
 520 static void analyse_disre(const char *voutfn,    int nframes,
 521                           real violaver[], real bounds[], int index[],
 522                           int pair[],      int nbounds,
 523                           const output_env_t oenv)
 524 {
 525     FILE   *vout;
 526     double  sum, sumt, sumaver;
 527     int     i, j;
 528
 529     /* Subtract bounds from distances, to calculate violations */
 530     calc_violations(violaver, violaver,
 531                     nbounds, pair, bounds, NULL, &sumt, &sumaver);
 532
 533 #ifdef DEBUG
 534     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 535             sumaver);
 536     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n",
 537             sumt);
 538 #endif
 539     vout = xvgropen(voutfn, "r\\S-3\\N average violations", "DR Index", "nm",
 540                     oenv);
 541     sum  = 0.0;
 542     sumt = 0.0;
 543     for (i = 0; (i < nbounds); i++)
 544     {
 545         /* Do ensemble averaging */
 546         sumaver = 0;
 547         for (j = pair[i]; (j < pair[i+1]); j++)
 548         {
 549             sumaver += sqr(violaver[j]/nframes);
 550         }
 551         sumaver = max(0.0, mypow(sumaver, minsixth)-bounds[i]);
 552
 553         sumt   += sumaver;
 554         sum     = max(sum, sumaver);
 555         fprintf(vout, "%10d  %10.5e\n", index[i], sumaver);
 556     }
 557 #ifdef DEBUG
 558     for (j = 0; (j < dr.ndr); j++)
 559     {
 560         fprintf(vout, "%10d  %10.5e\n", j, mypow(violaver[j]/nframes, minthird));
 561     }
 562 #endif
 563     gmx_ffclose(vout);
 564
 565     fprintf(stdout, "\nSum of violations averaged over simulation: %g nm\n",
 566             sumt);
 567     fprintf(stdout, "Largest violation averaged over simulation: %g nm\n\n", sum);
 568
 569     do_view(oenv, voutfn, "-graphtype bar");
 570 }
 571
 572 static void einstein_visco(const char *fn, const char *fni, int nsets,
 573                            int 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     gmx_ffclose(fp0);
 626     gmx_ffclose(fp1);
 627 }
 628
 629 typedef struct {
 630     gmx_int64_t     np;
 631     double          sum;
 632     double          sav;
 633     double          sav2;
 634 } ee_sum_t;
 635
 636 typedef struct {
 637     int             b;
 638     ee_sum_t        sum;
 639     gmx_int64_t     nst;
 640     gmx_int64_t     nst_min;
 641 } ener_ee_t;
 642
 643 static void clear_ee_sum(ee_sum_t *ees)
 644 {
 645     ees->sav  = 0;
 646     ees->sav2 = 0;
 647     ees->np   = 0;
 648     ees->sum  = 0;
 649 }
 650
 651 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
 652 {
 653     ees->np  += np;
 654     ees->sum += sum;
 655 }
 656
 657 static void add_ee_av(ee_sum_t *ees)
 658 {
 659     double av;
 660
 661     av         = ees->sum/ees->np;
 662     ees->sav  += av;
 663     ees->sav2 += av*av;
 664     ees->np    = 0;
 665     ees->sum   = 0;
 666 }
 667
 668 static double calc_ee2(int nb, ee_sum_t *ees)
 669 {
 670     return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
 671 }
 672
 673 static void set_ee_av(ener_ee_t *eee)
 674 {
 675     if (debug)
 676     {
 677         char buf[STEPSTRSIZE];
 678         fprintf(debug, "Storing average for err.est.: %s steps\n",
 679                 gmx_step_str(eee->nst, buf));
 680     }
 681     add_ee_av(&eee->sum);
 682     eee->b++;
 683     if (eee->b == 1 || eee->nst < eee->nst_min)
 684     {
 685         eee->nst_min = eee->nst;
 686     }
 687     eee->nst = 0;
 688 }
 689
 690 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
 691 {
 692     int             nb, i, f, nee;
 693     double          sum, sum2, sump, see2;
 694     gmx_int64_t     steps, np, p, bound_nb;
 695     enerdat_t      *ed;
 696     exactsum_t     *es;
 697     gmx_bool        bAllZero;
 698     double          x, sx, sy, sxx, sxy;
 699     ener_ee_t      *eee;
 700
 701     /* Check if we have exact statistics over all points */
 702     for (i = 0; i < nset; i++)
 703     {
 704         ed             = &edat->s[i];
 705         ed->bExactStat = FALSE;
 706         if (edat->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 nmol,
 972                                    char **leg, enerdata_t *edat,
 973                                    int nbmin, int nbmax)
 974 {
 975     int    i, j;
 976     double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, et, varet;
 977     double NANO3;
 978     enum {
 979         eVol, eEnth, eTemp, eEtot, eNR
 980     };
 981     const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
 982     int         ii[eNR];
 983
 984     NANO3 = NANO*NANO*NANO;
 985     if (!bDriftCorr)
 986     {
 987         fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
 988                 "for spurious drift in the graphs. Note that this is not\n"
 989                 "a substitute for proper equilibration and sampling!\n");
 990     }
 991     else
 992     {
 993         remove_drift(nset, nbmin, nbmax, dt, edat);
 994     }
 995     for (i = 0; (i < eNR); i++)
 996     {
 997         for (ii[i] = 0; (ii[i] < nset &&
 998                          (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
 999         {
1000             ;
1001         }
1002 /*        if (ii[i] < nset)
1003             fprintf(fp,"Found %s data.\n",my_ener[i]);
1004  */ }
1005     /* Compute it all! */
1006     alpha = kappa = cp = dcp = cv = NOTSET;
1007
1008     /* Temperature */
1009     tt = NOTSET;
1010     if (ii[eTemp] < nset)
1011     {
1012         tt    = edat->s[ii[eTemp]].av;
1013     }
1014     /* Volume */
1015     vv = varv = NOTSET;
1016     if ((ii[eVol] < nset) && (ii[eTemp] < nset))
1017     {
1018         vv    = edat->s[ii[eVol]].av*NANO3;
1019         varv  = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
1020         kappa = (varv/vv)/(BOLTZMANN*tt);
1021     }
1022     /* Enthalpy */
1023     hh = varh = NOTSET;
1024     if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
1025     {
1026         hh    = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
1027         varh  = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
1028         cp    = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
1029     }
1030     /* Total energy */
1031     et = varet = NOTSET;
1032     if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
1033     {
1034         /* Only compute cv in constant volume runs, which we can test
1035            by checking whether the enthalpy was computed.
1036          */
1037         et    = edat->s[ii[eEtot]].av;
1038         varet = sqr(edat->s[ii[eEtot]].rmsd);
1039         cv    = KILO*((varet/nmol)/(BOLTZ*tt*tt));
1040     }
1041     /* Alpha, dcp */
1042     if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
1043     {
1044         double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
1045         vh_sum = v_sum = h_sum = 0;
1046         for (j = 0; (j < edat->nframes); j++)
1047         {
1048             v       = edat->s[ii[eVol]].ener[j]*NANO3;
1049             h       = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
1050             v_sum  += v;
1051             h_sum  += h;
1052             vh_sum += (v*h);
1053         }
1054         vh_aver = vh_sum / edat->nframes;
1055         v_aver  = v_sum  / edat->nframes;
1056         h_aver  = h_sum  / edat->nframes;
1057         alpha   = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
1058         dcp     = (v_aver*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
1059     }
1060
1061     if (tt != NOTSET)
1062     {
1063         if (nmol < 2)
1064         {
1065             fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
1066                     nmol);
1067         }
1068         fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
1069         fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
1070         fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
1071         fprintf(fp, "please use the g_dos program.\n\n");
1072         fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
1073                 "a block-averaging error analysis (not implemented in g_energy yet)\n");
1074
1075         if (debug != NULL)
1076         {
1077             if (varv != NOTSET)
1078             {
1079                 fprintf(fp, "varv  =  %10g (m^6)\n", varv*AVOGADRO/nmol);
1080             }
1081         }
1082         if (vv != NOTSET)
1083         {
1084             fprintf(fp, "Volume                                   = %10g m^3/mol\n",
1085                     vv*AVOGADRO/nmol);
1086         }
1087         if (varh != NOTSET)
1088         {
1089             fprintf(fp, "Enthalpy                                 = %10g kJ/mol\n",
1090                     hh*AVOGADRO/(KILO*nmol));
1091         }
1092         if (alpha != NOTSET)
1093         {
1094             fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
1095                     alpha);
1096         }
1097         if (kappa != NOTSET)
1098         {
1099             fprintf(fp, "Isothermal Compressibility Kappa         = %10g (J/m^3)\n",
1100                     kappa);
1101             fprintf(fp, "Adiabatic bulk modulus                   = %10g (m^3/J)\n",
1102                     1.0/kappa);
1103         }
1104         if (cp != NOTSET)
1105         {
1106             fprintf(fp, "Heat capacity at constant pressure Cp    = %10g J/mol K\n",
1107                     cp);
1108         }
1109         if (cv != NOTSET)
1110         {
1111             fprintf(fp, "Heat capacity at constant volume Cv      = %10g J/mol K\n",
1112                     cv);
1113         }
1114         if (dcp != NOTSET)
1115         {
1116             fprintf(fp, "Cp-Cv                                    =  %10g J/mol K\n",
1117                     dcp);
1118         }
1119         please_cite(fp, "Allen1987a");
1120     }
1121     else
1122     {
1123         fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
1124     }
1125 }
1126
1127 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
1128                          gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
1129                          gmx_bool bVisco, const char *visfn, int nmol,
1130                          gmx_int64_t start_step, double start_t,
1131                          gmx_int64_t step, double t,
1132                          real reftemp,
1133                          enerdata_t *edat,
1134                          int nset, int set[], gmx_bool *bIsEner,
1135                          char **leg, gmx_enxnm_t *enm,
1136                          real Vaver, real ezero,
1137                          int nbmin, int nbmax,
1138                          const output_env_t oenv)
1139 {
1140     FILE           *fp;
1141     /* Check out the printed manual for equations! */
1142     double          Dt, aver, stddev, errest, delta_t, totaldrift;
1143     enerdata_t     *esum = NULL;
1144     real            xxx, integral, intBulk, Temp = 0, Pres = 0;
1145     real            sfrac, oldfrac, diffsum, diffav, fstep, pr_aver, pr_stddev, pr_errest;
1146     double          beta = 0, expE, expEtot, *fee = NULL;
1147     gmx_int64_t     nsteps;
1148     int             nexact, nnotexact;
1149     double          x1m, x1mk;
1150     int             i, j, k, nout;
1151     real            chi2;
1152     char            buf[256], eebuf[100];
1153
1154     nsteps  = step - start_step + 1;
1155     if (nsteps < 1)
1156     {
1157         fprintf(stdout, "Not enough steps (%s) for statistics\n",
1158                 gmx_step_str(nsteps, buf));
1159     }
1160     else
1161     {
1162         /* Calculate the time difference */
1163         delta_t = t - start_t;
1164
1165         fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1166                 gmx_step_str(nsteps, buf), start_t, t, nset);
1167
1168         calc_averages(nset, edat, nbmin, nbmax);
1169
1170         if (bSum)
1171         {
1172             esum = calc_sum(nset, edat, nbmin, nbmax);
1173         }
1174
1175         if (edat->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);
1399
1400             /* Now for bulk viscosity */
1401             strcpy(buf, "Bulk Viscosity");
1402             low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1403                             (edat->nframes+1)/2, &(eneset[11]), Dt,
1404                             eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1405
1406             factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1407             fp     = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1408             xvgr_legend(fp, asize(leg), leg, oenv);
1409
1410             /* Use trapezium rule for integration */
1411             integral = 0;
1412             intBulk  = 0;
1413             nout     = get_acfnout();
1414             if ((nout < 2) || (nout >= edat->nframes/2))
1415             {
1416                 nout = edat->nframes/2;
1417             }
1418             for (i = 1; (i < nout); i++)
1419             {
1420                 integral += 0.5*(eneset[0][i-1]  + eneset[0][i])*factor;
1421                 intBulk  += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1422                 fprintf(fp, "%10g  %10g  %10g\n", (i*Dt), integral, intBulk);
1423             }
1424             gmx_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         gmx_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_int64_t sum = 0;
1714         /* histograms */
1715         for (i = 0; i < fr->nblock; i++)
1716         {
1717             t_enxblock *blk = &(fr->block[i]);
1718             if (blk->id == enxDHHIST)
1719             {
1720                 double          foreign_lambda, dx;
1721                 gmx_int64_t     x0;
1722                 int             nhist, derivative;
1723
1724                 /* check the block types etc. */
1725                 if ( (blk->nsub < 2) ||
1726                      (blk->sub[0].type != xdr_datatype_double) ||
1727                      (blk->sub[1].type != xdr_datatype_int64) ||
1728                      (blk->sub[0].nr < 2)  ||
1729                      (blk->sub[1].nr < 2) )
1730                 {
1731                     gmx_fatal(FARGS, "Unexpected block data in file");
1732                 }
1733                 foreign_lambda = blk->sub[0].dval[0];
1734                 dx             = blk->sub[0].dval[1];
1735                 nhist          = blk->sub[1].lval[0];
1736                 derivative     = blk->sub[1].lval[1];
1737                 for (j = 0; j < nhist; j++)
1738                 {
1739                     const char *lg[1];
1740                     x0 = blk->sub[1].lval[2+j];
1741
1742                     if (!derivative)
1743                     {
1744                         sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1745                                 deltag, lambda, foreign_lambda,
1746                                 lambda, start_lambda);
1747                     }
1748                     else
1749                     {
1750                         sprintf(legend, "N(%s | %s=%g)",
1751                                 dhdl, lambda, start_lambda);
1752                     }
1753
1754                     lg[0] = legend;
1755                     xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1756                     setnr++;
1757                     for (k = 0; k < blk->sub[j+2].nr; k++)
1758                     {
1759                         int    hist;
1760                         double xmin, xmax;
1761
1762                         hist = blk->sub[j+2].ival[k];
1763                         xmin = (x0+k)*dx;
1764                         xmax = (x0+k+1)*dx;
1765                         fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1766                                 xmax, hist);
1767                         sum += hist;
1768                     }
1769                     /* multiple histogram data blocks in one histogram
1770                        mean that the second one is the reverse of the first one:
1771                        for dhdl derivatives, it's important to know both the
1772                        maximum and minimum values */
1773                     dx = -dx;
1774                 }
1775             }
1776         }
1777         (*samples) += (int)(sum/nblock_hist);
1778     }
1779     else
1780     {
1781         /* raw dh */
1782         int    len      = 0;
1783         char **setnames = NULL;
1784         int    nnames   = nblock_dh;
1785
1786         for (i = 0; i < fr->nblock; i++)
1787         {
1788             t_enxblock *blk = &(fr->block[i]);
1789             if (blk->id == enxDH)
1790             {
1791                 if (len == 0)
1792                 {
1793                     len = blk->sub[2].nr;
1794                 }
1795                 else
1796                 {
1797                     if (len != blk->sub[2].nr)
1798                     {
1799                         gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1800                     }
1801                 }
1802             }
1803         }
1804         (*samples) += len;
1805
1806         for (i = 0; i < len; i++)
1807         {
1808             double time = start_time + delta_time*i;
1809
1810             fprintf(*fp_dhdl, "%.4f ", time);
1811
1812             for (j = 0; j < fr->nblock; j++)
1813             {
1814                 t_enxblock *blk = &(fr->block[j]);
1815                 if (blk->id == enxDH)
1816                 {
1817                     double value;
1818                     if (blk->sub[2].type == xdr_datatype_float)
1819                     {
1820                         value = blk->sub[2].fval[i];
1821                     }
1822                     else
1823                     {
1824                         value = blk->sub[2].dval[i];
1825                     }
1826                     /* we need to decide which data type it is based on the count*/
1827
1828                     if (j == 1 && ir->bExpanded)
1829                     {
1830                         fprintf(*fp_dhdl, "%4d", (int)value);   /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
1831                     }
1832                     else
1833                     {
1834                         if (bDp)
1835                         {
1836                             fprintf(*fp_dhdl, " %#.12g", value);   /* print normal precision */
1837                         }
1838                         else
1839                         {
1840                             fprintf(*fp_dhdl, " %#.8g", value);   /* print normal precision */
1841                         }
1842                     }
1843                 }
1844             }
1845             fprintf(*fp_dhdl, "\n");
1846         }
1847     }
1848 }
1849
1850
1851 int gmx_energy(int argc, char *argv[])
1852 {
1853     const char        *desc[] = {
1854         "[THISMODULE] extracts energy components or distance restraint",
1855         "data from an energy file. The user is prompted to interactively",
1856         "select the desired energy terms.[PAR]",
1857
1858         "Average, RMSD, and drift are calculated with full precision from the",
1859         "simulation (see printed manual). Drift is calculated by performing",
1860         "a least-squares fit of the data to a straight line. The reported total drift",
1861         "is the difference of the fit at the first and last point.",
1862         "An error estimate of the average is given based on a block averages",
1863         "over 5 blocks using the full-precision averages. The error estimate",
1864         "can be performed over multiple block lengths with the options",
1865         "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1866         "[BB]Note[bb] that in most cases the energy files contains averages over all",
1867         "MD steps, or over many more points than the number of frames in",
1868         "energy file. This makes the [THISMODULE] statistics output more accurate",
1869         "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1870         "file, the statistics mentioned above are simply over the single, per-frame",
1871         "energy values.[PAR]",
1872
1873         "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1874
1875         "Some fluctuation-dependent properties can be calculated provided",
1876         "the correct energy terms are selected, and that the command line option",
1877         "[TT]-fluct_props[tt] is given. The following properties",
1878         "will be computed:[BR]",
1879         "Property                        Energy terms needed[BR]",
1880         "---------------------------------------------------[BR]",
1881         "Heat capacity C[SUB]p[sub] (NPT sims):    Enthalpy, Temp     [BR]",
1882         "Heat capacity C[SUB]v[sub] (NVT sims):    Etot, Temp         [BR]",
1883         "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1884         "Isothermal compressibility:     Vol, Temp          [BR]",
1885         "Adiabatic bulk modulus:         Vol, Temp          [BR]",
1886         "---------------------------------------------------[BR]",
1887         "You always need to set the number of molecules [TT]-nmol[tt].",
1888         "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1889         "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]"
1890         "When the [TT]-viol[tt] option is set, the time averaged",
1891         "violations are plotted and the running time-averaged and",
1892         "instantaneous sum of violations are recalculated. Additionally",
1893         "running time-averaged and instantaneous distances between",
1894         "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1895
1896         "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1897         "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1898         "The first two options plot the orientation, the last three the",
1899         "deviations of the orientations from the experimental values.",
1900         "The options that end on an 'a' plot the average over time",
1901         "as a function of restraint. The options that end on a 't'",
1902         "prompt the user for restraint label numbers and plot the data",
1903         "as a function of time. Option [TT]-odr[tt] plots the RMS",
1904         "deviation as a function of restraint.",
1905         "When the run used time or ensemble averaged orientation restraints,",
1906         "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1907         "not ensemble-averaged orientations and deviations instead of",
1908         "the time and ensemble averages.[PAR]",
1909
1910         "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1911         "tensor for each orientation restraint experiment. With option",
1912         "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1913
1914         "Option [TT]-odh[tt] extracts and plots the free energy data",
1915         "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1916         "from the [TT]ener.edr[tt] file.[PAR]",
1917
1918         "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1919         "difference with an ideal gas state: [BR]",
1920         "  [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]",
1921         "  [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]",
1922         "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1923         "the average is over the ensemble (or time in a trajectory).",
1924         "Note that this is in principle",
1925         "only correct when averaging over the whole (Boltzmann) ensemble",
1926         "and using the potential energy. This also allows for an entropy",
1927         "estimate using:[BR]",
1928         "  [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]",
1929         "  [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",
1930         "[PAR]",
1931
1932         "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1933         "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1934         "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1935         "files, and the average is over the ensemble A. The running average",
1936         "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1937         "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1938
1939     };
1940     static gmx_bool    bSum    = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1941     static gmx_bool    bDp     = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1942     static int         skip    = 0, nmol = 1, nbmin = 5, nbmax = 5;
1943     static real        reftemp = 300.0, ezero = 0;
1944     t_pargs            pa[]    = {
1945         { "-fee",   FALSE, etBOOL,  {&bFee},
1946           "Do a free energy estimate" },
1947         { "-fetemp", FALSE, etREAL, {&reftemp},
1948           "Reference temperature for free energy calculation" },
1949         { "-zero", FALSE, etREAL, {&ezero},
1950           "Subtract a zero-point energy" },
1951         { "-sum",  FALSE, etBOOL, {&bSum},
1952           "Sum the energy terms selected rather than display them all" },
1953         { "-dp",   FALSE, etBOOL, {&bDp},
1954           "Print energies in high precision" },
1955         { "-nbmin", FALSE, etINT, {&nbmin},
1956           "Minimum number of blocks for error estimate" },
1957         { "-nbmax", FALSE, etINT, {&nbmax},
1958           "Maximum number of blocks for error estimate" },
1959         { "-mutot", FALSE, etBOOL, {&bMutot},
1960           "Compute the total dipole moment from the components" },
1961         { "-skip", FALSE, etINT,  {&skip},
1962           "Skip number of frames between data points" },
1963         { "-aver", FALSE, etBOOL, {&bPrAll},
1964           "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1965         { "-nmol", FALSE, etINT,  {&nmol},
1966           "Number of molecules in your sample: the energies are divided by this number" },
1967         { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1968           "Compute properties based on energy fluctuations, like heat capacity" },
1969         { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1970           "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1971         { "-fluc", FALSE, etBOOL, {&bFluct},
1972           "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1973         { "-orinst", FALSE, etBOOL, {&bOrinst},
1974           "Analyse instantaneous orientation data" },
1975         { "-ovec", FALSE, etBOOL, {&bOvec},
1976           "Also plot the eigenvectors with [TT]-oten[tt]" }
1977     };
1978     const char       * drleg[] = {
1979         "Running average",
1980         "Instantaneous"
1981     };
1982     static const char *setnm[] = {
1983         "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1984         "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1985         "Volume",  "Pressure"
1986     };
1987
1988     FILE              *out     = NULL, *fp_pairs = NULL, *fort = NULL, *fodt = NULL, *foten = NULL;
1989     FILE              *fp_dhdl = NULL;
1990     FILE             **drout;
1991     ener_file_t        fp;
1992     int                timecheck = 0;
1993     gmx_mtop_t         mtop;
1994     gmx_localtop_t    *top = NULL;
1995     t_inputrec         ir;
1996     t_energy         **ee;
1997     enerdata_t         edat;
1998     gmx_enxnm_t       *enm = NULL;
1999     t_enxframe        *frame, *fr = NULL;
2000     int                cur = 0;
2001 #define NEXT (1-cur)
2002     int                nre, teller, teller_disre, nfr;
2003     gmx_int64_t        start_step;
2004     int                nor = 0, nex = 0, norfr = 0, enx_i = 0;
2005     real               start_t;
2006     real              *bounds  = NULL, *violaver = NULL, *oobs = NULL, *orient = NULL, *odrms = NULL;
2007     int               *index   = NULL, *pair = NULL, norsel = 0, *orsel = NULL, *or_label = NULL;
2008     int                nbounds = 0, npairs;
2009     gmx_bool           bDisRe, bDRAll, bORA, bORT, bODA, bODR, bODT, bORIRE, bOTEN, bDHDL;
2010     gmx_bool           bFoundStart, bCont, bEDR, bVisco;
2011     double             sum, sumaver, sumt, ener, dbl;
2012     double            *time = NULL;
2013     real               Vaver;
2014     int               *set     = NULL, i, j, k, nset, sss;
2015     gmx_bool          *bIsEner = NULL;
2016     char             **pairleg, **odtleg, **otenleg;
2017     char             **leg = NULL;
2018     char             **nms;
2019     char              *anm_j, *anm_k, *resnm_j, *resnm_k;
2020     int                resnr_j, resnr_k;
2021     const char        *orinst_sub = "@ subtitle \"instantaneous\"\n";
2022     char               buf[256];
2023     output_env_t       oenv;
2024     t_enxblock        *blk       = NULL;
2025     t_enxblock        *blk_disre = NULL;
2026     int                ndisre    = 0;
2027     int                dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
2028
2029     t_filenm           fnm[] = {
2030         { efEDR, "-f",    NULL,      ffREAD  },
2031         { efEDR, "-f2",   NULL,      ffOPTRD },
2032         { efTPX, "-s",    NULL,      ffOPTRD },
2033         { efXVG, "-o",    "energy",  ffWRITE },
2034         { efXVG, "-viol", "violaver", ffOPTWR },
2035         { efXVG, "-pairs", "pairs",   ffOPTWR },
2036         { efXVG, "-ora",  "orienta", ffOPTWR },
2037         { efXVG, "-ort",  "orientt", ffOPTWR },
2038         { efXVG, "-oda",  "orideva", ffOPTWR },
2039         { efXVG, "-odr",  "oridevr", ffOPTWR },
2040         { efXVG, "-odt",  "oridevt", ffOPTWR },
2041         { efXVG, "-oten", "oriten",  ffOPTWR },
2042         { efXVG, "-corr", "enecorr", ffOPTWR },
2043         { efXVG, "-vis",  "visco",   ffOPTWR },
2044         { efXVG, "-ravg", "runavgdf", ffOPTWR },
2045         { efXVG, "-odh",  "dhdl", ffOPTWR }
2046     };
2047 #define NFILE asize(fnm)
2048     int                npargs;
2049     t_pargs           *ppa;
2050
2051     npargs = asize(pa);
2052     ppa    = add_acf_pargs(&npargs, pa);
2053     if (!parse_common_args(&argc, argv,
2054                            PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END,
2055                            NFILE, fnm, npargs, ppa, asize(desc), desc, 0, NULL, &oenv))
2056     {
2057         return 0;
2058     }
2059
2060     bDRAll = opt2bSet("-pairs", NFILE, fnm);
2061     bDisRe = opt2bSet("-viol", NFILE, fnm) || bDRAll;
2062     bORA   = opt2bSet("-ora", NFILE, fnm);
2063     bORT   = opt2bSet("-ort", NFILE, fnm);
2064     bODA   = opt2bSet("-oda", NFILE, fnm);
2065     bODR   = opt2bSet("-odr", NFILE, fnm);
2066     bODT   = opt2bSet("-odt", NFILE, fnm);
2067     bORIRE = bORA || bORT || bODA || bODR || bODT;
2068     bOTEN  = opt2bSet("-oten", NFILE, fnm);
2069     bDHDL  = opt2bSet("-odh", NFILE, fnm);
2070
2071     nset = 0;
2072
2073     snew(frame, 2);
2074     fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
2075     do_enxnms(fp, &nre, &enm);
2076
2077     Vaver = -1;
2078
2079     bVisco = opt2bSet("-vis", NFILE, fnm);
2080
2081     if ((!bDisRe) && (!bDHDL))
2082     {
2083         if (bVisco)
2084         {
2085             nset = asize(setnm);
2086             snew(set, nset);
2087             /* This is nasty code... To extract Pres tensor, Volume and Temperature */
2088             for (j = 0; j < nset; j++)
2089             {
2090                 for (i = 0; i < nre; i++)
2091                 {
2092                     if (strstr(enm[i].name, setnm[j]))
2093                     {
2094                         set[j] = i;
2095                         break;
2096                     }
2097                 }
2098                 if (i == nre)
2099                 {
2100                     if (gmx_strcasecmp(setnm[j], "Volume") == 0)
2101                     {
2102                         printf("Enter the box volume (" unit_volume "): ");
2103                         if (1 != scanf("%lf", &dbl))
2104                         {
2105                             gmx_fatal(FARGS, "Error reading user input");
2106                         }
2107                         Vaver = dbl;
2108                     }
2109                     else
2110                     {
2111                         gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
2112                                   setnm[j]);
2113                     }
2114                 }
2115             }
2116         }
2117         else
2118         {
2119             set = select_by_name(nre, enm, &nset);
2120         }
2121         /* Print all the different units once */
2122         sprintf(buf, "(%s)", enm[set[0]].unit);
2123         for (i = 1; i < nset; i++)
2124         {
2125             for (j = 0; j < i; j++)
2126             {
2127                 if (strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
2128                 {
2129                     break;
2130                 }
2131             }
2132             if (j == i)
2133             {
2134                 strcat(buf, ", (");
2135                 strcat(buf, enm[set[i]].unit);
2136                 strcat(buf, ")");
2137             }
2138         }
2139         out = xvgropen(opt2fn("-o", NFILE, fnm), "Gromacs Energies", "Time (ps)", buf,
2140                        oenv);
2141
2142         snew(leg, nset+1);
2143         for (i = 0; (i < nset); i++)
2144         {
2145             leg[i] = enm[set[i]].name;
2146         }
2147         if (bSum)
2148         {
2149             leg[nset] = gmx_strdup("Sum");
2150             xvgr_legend(out, nset+1, (const char**)leg, oenv);
2151         }
2152         else
2153         {
2154             xvgr_legend(out, nset, (const char**)leg, oenv);
2155         }
2156
2157         snew(bIsEner, nset);
2158         for (i = 0; (i < nset); i++)
2159         {
2160             bIsEner[i] = FALSE;
2161             for (j = 0; (j <= F_ETOT); j++)
2162             {
2163                 bIsEner[i] = bIsEner[i] ||
2164                     (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
2165             }
2166         }
2167
2168         if (bPrAll && nset > 1)
2169         {
2170             gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
2171         }
2172
2173         time = NULL;
2174
2175         if (bORIRE || bOTEN)
2176         {
2177             get_orires_parms(ftp2fn(efTPX, NFILE, fnm), &nor, &nex, &or_label, &oobs);
2178         }
2179
2180         if (bORIRE)
2181         {
2182             if (bOrinst)
2183             {
2184                 enx_i = enxORI;
2185             }
2186             else
2187             {
2188                 enx_i = enxOR;
2189             }
2190
2191             if (bORA || bODA)
2192             {
2193                 snew(orient, nor);
2194             }
2195             if (bODR)
2196             {
2197                 snew(odrms, nor);
2198             }
2199             if (bORT || bODT)
2200             {
2201                 fprintf(stderr, "Select the orientation restraint labels you want (-1 is all)\n");
2202                 fprintf(stderr, "End your selection with 0\n");
2203                 j     = -1;
2204                 orsel = NULL;
2205                 do
2206                 {
2207                     j++;
2208                     srenew(orsel, j+1);
2209                     if (1 != scanf("%d", &(orsel[j])))
2210                     {
2211                         gmx_fatal(FARGS, "Error reading user input");
2212                     }
2213                 }
2214                 while (orsel[j] > 0);
2215                 if (orsel[0] == -1)
2216                 {
2217                     fprintf(stderr, "Selecting all %d orientation restraints\n", nor);
2218                     norsel = nor;
2219                     srenew(orsel, nor);
2220                     for (i = 0; i < nor; i++)
2221                     {
2222                         orsel[i] = i;
2223                     }
2224                 }
2225                 else
2226                 {
2227                     /* Build the selection */
2228                     norsel = 0;
2229                     for (i = 0; i < j; i++)
2230                     {
2231                         for (k = 0; k < nor; k++)
2232                         {
2233                             if (or_label[k] == orsel[i])
2234                             {
2235                                 orsel[norsel] = k;
2236                                 norsel++;
2237                                 break;
2238                             }
2239                         }
2240                         if (k == nor)
2241                         {
2242                             fprintf(stderr, "Orientation restraint label %d not found\n",
2243                                     orsel[i]);
2244                         }
2245                     }
2246                 }
2247                 snew(odtleg, norsel);
2248                 for (i = 0; i < norsel; i++)
2249                 {
2250                     snew(odtleg[i], 256);
2251                     sprintf(odtleg[i], "%d", or_label[orsel[i]]);
2252                 }
2253                 if (bORT)
2254                 {
2255                     fort = xvgropen(opt2fn("-ort", NFILE, fnm), "Calculated orientations",
2256                                     "Time (ps)", "", oenv);
2257                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2258                     {
2259                         fprintf(fort, "%s", orinst_sub);
2260                     }
2261                     xvgr_legend(fort, norsel, (const char**)odtleg, oenv);
2262                 }
2263                 if (bODT)
2264                 {
2265                     fodt = xvgropen(opt2fn("-odt", NFILE, fnm),
2266                                     "Orientation restraint deviation",
2267                                     "Time (ps)", "", oenv);
2268                     if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2269                     {
2270                         fprintf(fodt, "%s", orinst_sub);
2271                     }
2272                     xvgr_legend(fodt, norsel, (const char**)odtleg, oenv);
2273                 }
2274             }
2275         }
2276         if (bOTEN)
2277         {
2278             foten = xvgropen(opt2fn("-oten", NFILE, fnm),
2279                              "Order tensor", "Time (ps)", "", oenv);
2280             snew(otenleg, bOvec ? nex*12 : nex*3);
2281             for (i = 0; i < nex; i++)
2282             {
2283                 for (j = 0; j < 3; j++)
2284                 {
2285                     sprintf(buf, "eig%d", j+1);
2286                     otenleg[(bOvec ? 12 : 3)*i+j] = gmx_strdup(buf);
2287                 }
2288                 if (bOvec)
2289                 {
2290                     for (j = 0; j < 9; j++)
2291                     {
2292                         sprintf(buf, "vec%d%s", j/3+1, j%3 == 0 ? "x" : (j%3 == 1 ? "y" : "z"));
2293                         otenleg[12*i+3+j] = gmx_strdup(buf);
2294                     }
2295                 }
2296             }
2297             xvgr_legend(foten, bOvec ? nex*12 : nex*3, (const char**)otenleg, oenv);
2298         }
2299     }
2300     else if (bDisRe)
2301     {
2302         nbounds = get_bounds(ftp2fn(efTPX, NFILE, fnm), &bounds, &index, &pair, &npairs,
2303                              &mtop, &top, &ir);
2304         snew(violaver, npairs);
2305         out = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations",
2306                        "Time (ps)", "nm", oenv);
2307         xvgr_legend(out, 2, drleg, oenv);
2308         if (bDRAll)
2309         {
2310             fp_pairs = xvgropen(opt2fn("-pairs", NFILE, fnm), "Pair Distances",
2311                                 "Time (ps)", "Distance (nm)", oenv);
2312             if (output_env_get_print_xvgr_codes(oenv))
2313             {
2314                 fprintf(fp_pairs, "@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2315                         ir.dr_tau);
2316             }
2317         }
2318     }
2319     else if (bDHDL)
2320     {
2321         get_dhdl_parms(ftp2fn(efTPX, NFILE, fnm), &ir);
2322     }
2323
2324     /* Initiate energies and set them to zero */
2325     edat.nsteps  = 0;
2326     edat.npoints = 0;
2327     edat.nframes = 0;
2328     edat.step    = NULL;
2329     edat.steps   = NULL;
2330     edat.points  = NULL;
2331     snew(edat.s, nset);
2332
2333     /* Initiate counters */
2334     teller       = 0;
2335     teller_disre = 0;
2336     bFoundStart  = FALSE;
2337     start_step   = 0;
2338     start_t      = 0;
2339     do
2340     {
2341         /* This loop searches for the first frame (when -b option is given),
2342          * or when this has been found it reads just one energy frame
2343          */
2344         do
2345         {
2346             bCont = do_enx(fp, &(frame[NEXT]));
2347             if (bCont)
2348             {
2349                 timecheck = check_times(frame[NEXT].t);
2350             }
2351         }
2352         while (bCont && (timecheck < 0));
2353
2354         if ((timecheck == 0) && bCont)
2355         {
2356             /* We read a valid frame, so we can use it */
2357             fr = &(frame[NEXT]);
2358
2359             if (fr->nre > 0)
2360             {
2361                 /* The frame contains energies, so update cur */
2362                 cur  = NEXT;
2363
2364                 if (edat.nframes % 1000 == 0)
2365                 {
2366                     srenew(edat.step, edat.nframes+1000);
2367                     memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
2368                     srenew(edat.steps, edat.nframes+1000);
2369                     memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
2370                     srenew(edat.points, edat.nframes+1000);
2371                     memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
2372
2373                     for (i = 0; i < nset; i++)
2374                     {
2375                         srenew(edat.s[i].ener, edat.nframes+1000);
2376                         memset(&(edat.s[i].ener[edat.nframes]), 0,
2377                                1000*sizeof(edat.s[i].ener[0]));
2378                         srenew(edat.s[i].es, edat.nframes+1000);
2379                         memset(&(edat.s[i].es[edat.nframes]), 0,
2380                                1000*sizeof(edat.s[i].es[0]));
2381                     }
2382                 }
2383
2384                 nfr            = edat.nframes;
2385                 edat.step[nfr] = fr->step;
2386
2387                 if (!bFoundStart)
2388                 {
2389                     bFoundStart = TRUE;
2390                     /* Initiate the previous step data */
2391                     start_step = fr->step;
2392                     start_t    = fr->t;
2393                     /* Initiate the energy sums */
2394                     edat.steps[nfr]  = 1;
2395                     edat.points[nfr] = 1;
2396                     for (i = 0; i < nset; i++)
2397                     {
2398                         sss                    = set[i];
2399                         edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2400                         edat.s[i].es[nfr].sum2 = 0;
2401                     }
2402                     edat.nsteps  = 1;
2403                     edat.npoints = 1;
2404                 }
2405                 else
2406                 {
2407                     edat.steps[nfr] = fr->nsteps;
2408                     {
2409                         if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2410                         {
2411                             if (fr->nsum <= 1)
2412                             {
2413                                 edat.points[nfr] = 1;
2414                                 for (i = 0; i < nset; i++)
2415                                 {
2416                                     sss                    = set[i];
2417                                     edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2418                                     edat.s[i].es[nfr].sum2 = 0;
2419                                 }
2420                                 edat.npoints += 1;
2421                             }
2422                             else
2423                             {
2424                                 edat.points[nfr] = fr->nsum;
2425                                 for (i = 0; i < nset; i++)
2426                                 {
2427                                     sss                    = set[i];
2428                                     edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2429                                     edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2430                                 }
2431                                 edat.npoints += fr->nsum;
2432                             }
2433                         }
2434                         else
2435                         {
2436                             /* The interval does not match fr->nsteps:
2437                              * can not do exact averages.
2438                              */
2439                             edat.npoints = 0;
2440                         }
2441                         edat.nsteps = fr->step - start_step + 1;
2442                     }
2443                 }
2444                 for (i = 0; i < nset; i++)
2445                 {
2446                     edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2447                 }
2448             }
2449             /*
2450              * Define distance restraint legends. Can only be done after
2451              * the first frame has been read... (Then we know how many there are)
2452              */
2453             blk_disre = find_block_id_enxframe(fr, enxDISRE, NULL);
2454             if (bDisRe && bDRAll && !leg && blk_disre)
2455             {
2456                 t_iatom   *fa;
2457                 t_iparams *ip;
2458
2459                 fa = top->idef.il[F_DISRES].iatoms;
2460                 ip = top->idef.iparams;
2461                 if (blk_disre->nsub != 2 ||
2462                     (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2463                 {
2464                     gmx_incons("Number of disre sub-blocks not equal to 2");
2465                 }
2466
2467                 ndisre = blk_disre->sub[0].nr;
2468                 if (ndisre != top->idef.il[F_DISRES].nr/3)
2469                 {
2470                     gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2471                               ndisre, top->idef.il[F_DISRES].nr/3);
2472                 }
2473                 snew(pairleg, ndisre);
2474                 for (i = 0; i < ndisre; i++)
2475                 {
2476                     snew(pairleg[i], 30);
2477                     j = fa[3*i+1];
2478                     k = fa[3*i+2];
2479                     gmx_mtop_atominfo_global(&mtop, j, &anm_j, &resnr_j, &resnm_j);
2480                     gmx_mtop_atominfo_global(&mtop, k, &anm_k, &resnr_k, &resnm_k);
2481                     sprintf(pairleg[i], "%d %s %d %s (%d)",
2482                             resnr_j, anm_j, resnr_k, anm_k,
2483                             ip[fa[3*i]].disres.label);
2484                 }
2485                 set = select_it(ndisre, pairleg, &nset);
2486                 snew(leg, 2*nset);
2487                 for (i = 0; (i < nset); i++)
2488                 {
2489                     snew(leg[2*i], 32);
2490                     sprintf(leg[2*i],  "a %s", pairleg[set[i]]);
2491                     snew(leg[2*i+1], 32);
2492                     sprintf(leg[2*i+1], "i %s", pairleg[set[i]]);
2493                 }
2494                 xvgr_legend(fp_pairs, 2*nset, (const char**)leg, oenv);
2495             }
2496
2497             /*
2498              * Store energies for analysis afterwards...
2499              */
2500             if (!bDisRe && !bDHDL && (fr->nre > 0))
2501             {
2502                 if (edat.nframes % 1000 == 0)
2503                 {
2504                     srenew(time, edat.nframes+1000);
2505                 }
2506                 time[edat.nframes] = fr->t;
2507                 edat.nframes++;
2508             }
2509             /*
2510              * Printing time, only when we do not want to skip frames
2511              */
2512             if (!skip || teller % skip == 0)
2513             {
2514                 if (bDisRe)
2515                 {
2516                     /*******************************************
2517                      * D I S T A N C E   R E S T R A I N T S
2518                      *******************************************/
2519                     if (ndisre > 0)
2520                     {
2521  #ifndef GMX_DOUBLE
2522                         float  *disre_rt     =     blk_disre->sub[0].fval;
2523                         float  *disre_rm3tav = blk_disre->sub[1].fval;
2524  #else
2525                         double *disre_rt     =     blk_disre->sub[0].dval;
2526                         double *disre_rm3tav = blk_disre->sub[1].dval;
2527  #endif
2528
2529                         print_time(out, fr->t);
2530                         if (violaver == NULL)
2531                         {
2532                             snew(violaver, ndisre);
2533                         }
2534
2535                         /* Subtract bounds from distances, to calculate violations */
2536                         calc_violations(disre_rt, disre_rm3tav,
2537                                         nbounds, pair, bounds, violaver, &sumt, &sumaver);
2538
2539                         fprintf(out, "  %8.4f  %8.4f\n", sumaver, sumt);
2540                         if (bDRAll)
2541                         {
2542                             print_time(fp_pairs, fr->t);
2543                             for (i = 0; (i < nset); i++)
2544                             {
2545                                 sss = set[i];
2546                                 fprintf(fp_pairs, "  %8.4f", mypow(disre_rm3tav[sss], minthird));
2547                                 fprintf(fp_pairs, "  %8.4f", disre_rt[sss]);
2548                             }
2549                             fprintf(fp_pairs, "\n");
2550                         }
2551                         teller_disre++;
2552                     }
2553                 }
2554                 else if (bDHDL)
2555                 {
2556                     do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2557                 }
2558
2559                 /*******************************************
2560                  * E N E R G I E S
2561                  *******************************************/
2562                 else
2563                 {
2564                     if (fr->nre > 0)
2565                     {
2566                         if (bPrAll)
2567                         {
2568                             /* We skip frames with single points (usually only the first frame),
2569                              * since they would result in an average plot with outliers.
2570                              */
2571                             if (fr->nsum > 1)
2572                             {
2573                                 print_time(out, fr->t);
2574                                 print1(out, bDp, fr->ener[set[0]].e);
2575                                 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2576                                 print1(out, bDp, sqrt(fr->ener[set[0]].eav/fr->nsum));
2577                                 fprintf(out, "\n");
2578                             }
2579                         }
2580                         else
2581                         {
2582                             print_time(out, fr->t);
2583                             if (bSum)
2584                             {
2585                                 sum = 0;
2586                                 for (i = 0; i < nset; i++)
2587                                 {
2588                                     sum += fr->ener[set[i]].e;
2589                                 }
2590                                 print1(out, bDp, sum/nmol-ezero);
2591                             }
2592                             else
2593                             {
2594                                 for (i = 0; (i < nset); i++)
2595                                 {
2596                                     if (bIsEner[i])
2597                                     {
2598                                         print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2599                                     }
2600                                     else
2601                                     {
2602                                         print1(out, bDp, fr->ener[set[i]].e);
2603                                     }
2604                                 }
2605                             }
2606                             fprintf(out, "\n");
2607                         }
2608                     }
2609                     blk = find_block_id_enxframe(fr, enx_i, NULL);
2610                     if (bORIRE && blk)
2611                     {
2612 #ifndef GMX_DOUBLE
2613                         xdr_datatype dt = xdr_datatype_float;
2614 #else
2615                         xdr_datatype dt = xdr_datatype_double;
2616 #endif
2617                         real        *vals;
2618
2619                         if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2620                         {
2621                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2622                         }
2623 #ifndef GMX_DOUBLE
2624                         vals = blk->sub[0].fval;
2625 #else
2626                         vals = blk->sub[0].dval;
2627 #endif
2628
2629                         if (blk->sub[0].nr != (size_t)nor)
2630                         {
2631                             gmx_fatal(FARGS, "Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2632                         }
2633                         if (bORA || bODA)
2634                         {
2635                             for (i = 0; i < nor; i++)
2636                             {
2637                                 orient[i] += vals[i];
2638                             }
2639                         }
2640                         if (bODR)
2641                         {
2642                             for (i = 0; i < nor; i++)
2643                             {
2644                                 odrms[i] += sqr(vals[i]-oobs[i]);
2645                             }
2646                         }
2647                         if (bORT)
2648                         {
2649                             fprintf(fort, "  %10f", fr->t);
2650                             for (i = 0; i < norsel; i++)
2651                             {
2652                                 fprintf(fort, " %g", vals[orsel[i]]);
2653                             }
2654                             fprintf(fort, "\n");
2655                         }
2656                         if (bODT)
2657                         {
2658                             fprintf(fodt, "  %10f", fr->t);
2659                             for (i = 0; i < norsel; i++)
2660                             {
2661                                 fprintf(fodt, " %g", vals[orsel[i]]-oobs[orsel[i]]);
2662                             }
2663                             fprintf(fodt, "\n");
2664                         }
2665                         norfr++;
2666                     }
2667                     blk = find_block_id_enxframe(fr, enxORT, NULL);
2668                     if (bOTEN && blk)
2669                     {
2670 #ifndef GMX_DOUBLE
2671                         xdr_datatype dt = xdr_datatype_float;
2672 #else
2673                         xdr_datatype dt = xdr_datatype_double;
2674 #endif
2675                         real        *vals;
2676
2677                         if ( (blk->nsub != 1) || (blk->sub[0].type != dt) )
2678                         {
2679                             gmx_fatal(FARGS, "Orientational restraints read in incorrectly");
2680                         }
2681 #ifndef GMX_DOUBLE
2682                         vals = blk->sub[0].fval;
2683 #else
2684                         vals = blk->sub[0].dval;
2685 #endif
2686
2687                         if (blk->sub[0].nr != (size_t)(nex*12))
2688                         {
2689                             gmx_fatal(FARGS, "Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2690                                       blk->sub[0].nr/12, nex);
2691                         }
2692                         fprintf(foten, "  %10f", fr->t);
2693                         for (i = 0; i < nex; i++)
2694                         {
2695                             for (j = 0; j < (bOvec ? 12 : 3); j++)
2696                             {
2697                                 fprintf(foten, " %g", vals[i*12+j]);
2698                             }
2699                         }
2700                         fprintf(foten, "\n");
2701                     }
2702                 }
2703             }
2704             teller++;
2705         }
2706     }
2707     while (bCont && (timecheck == 0));
2708
2709     fprintf(stderr, "\n");
2710     close_enx(fp);
2711     if (out)
2712     {
2713         gmx_ffclose(out);
2714     }
2715
2716     if (bDRAll)
2717     {
2718         gmx_ffclose(fp_pairs);
2719     }
2720
2721     if (bORT)
2722     {
2723         gmx_ffclose(fort);
2724     }
2725     if (bODT)
2726     {
2727         gmx_ffclose(fodt);
2728     }
2729     if (bORA)
2730     {
2731         out = xvgropen(opt2fn("-ora", NFILE, fnm),
2732                        "Average calculated orientations",
2733                        "Restraint label", "", oenv);
2734         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2735         {
2736             fprintf(out, "%s", orinst_sub);
2737         }
2738         for (i = 0; i < nor; i++)
2739         {
2740             fprintf(out, "%5d  %g\n", or_label[i], orient[i]/norfr);
2741         }
2742         gmx_ffclose(out);
2743     }
2744     if (bODA)
2745     {
2746         out = xvgropen(opt2fn("-oda", NFILE, fnm),
2747                        "Average restraint deviation",
2748                        "Restraint label", "", oenv);
2749         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2750         {
2751             fprintf(out, "%s", orinst_sub);
2752         }
2753         for (i = 0; i < nor; i++)
2754         {
2755             fprintf(out, "%5d  %g\n", or_label[i], orient[i]/norfr-oobs[i]);
2756         }
2757         gmx_ffclose(out);
2758     }
2759     if (bODR)
2760     {
2761         out = xvgropen(opt2fn("-odr", NFILE, fnm),
2762                        "RMS orientation restraint deviations",
2763                        "Restraint label", "", oenv);
2764         if (bOrinst && output_env_get_print_xvgr_codes(oenv))
2765         {
2766             fprintf(out, "%s", orinst_sub);
2767         }
2768         for (i = 0; i < nor; i++)
2769         {
2770             fprintf(out, "%5d  %g\n", or_label[i], sqrt(odrms[i]/norfr));
2771         }
2772         gmx_ffclose(out);
2773     }
2774     if (bOTEN)
2775     {
2776         gmx_ffclose(foten);
2777     }
2778
2779     if (bDisRe)
2780     {
2781         analyse_disre(opt2fn("-viol", NFILE, fnm),
2782                       teller_disre, violaver, bounds, index, pair, nbounds, oenv);
2783     }
2784     else if (bDHDL)
2785     {
2786         if (fp_dhdl)
2787         {
2788             gmx_ffclose(fp_dhdl);
2789             printf("\n\nWrote %d lambda values with %d samples as ",
2790                    dh_lambdas, dh_samples);
2791             if (dh_hists > 0)
2792             {
2793                 printf("%d dH histograms ", dh_hists);
2794             }
2795             if (dh_blocks > 0)
2796             {
2797                 printf("%d dH data blocks ", dh_blocks);
2798             }
2799             printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2800
2801         }
2802         else
2803         {
2804             gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2805         }
2806
2807     }
2808     else
2809     {
2810         double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2811         analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2812                      bFee, bSum, opt2parg_bSet("-nmol", npargs, ppa),
2813                      bVisco, opt2fn("-vis", NFILE, fnm),
2814                      nmol,
2815                      start_step, start_t, frame[cur].step, frame[cur].t,
2816                      reftemp, &edat,
2817                      nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2818                      oenv);
2819         if (bFluctProps)
2820         {
2821             calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2822                                    nbmin, nbmax);
2823         }
2824     }
2825     if (opt2bSet("-f2", NFILE, fnm))
2826     {
2827         fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2828             reftemp, nset, set, leg, &edat, time, oenv);
2829     }
2830
2831     {
2832         const char *nxy = "-nxy";
2833
2834         do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2835         do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2836         do_view(oenv, opt2fn_null("-ora", NFILE, fnm), nxy);
2837         do_view(oenv, opt2fn_null("-ort", NFILE, fnm), nxy);
2838         do_view(oenv, opt2fn_null("-oda", NFILE, fnm), nxy);
2839         do_view(oenv, opt2fn_null("-odr", NFILE, fnm), nxy);
2840         do_view(oenv, opt2fn_null("-odt", NFILE, fnm), nxy);
2841         do_view(oenv, opt2fn_null("-oten", NFILE, fnm), nxy);
2842         do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);
2843     }
2844
2845     return 0;
2846 }