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