src/gromacs/gmxana/gmx_dos.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2011,2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "gmxpre.h"
  36
  37 #include "config.h"
  38
  39 #include <math.h>
  40 #include <stdio.h>
  41 #include <stdlib.h>
  42 #include <string.h>
  43
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/legacyheaders/copyrite.h"
  46 #include "gromacs/utility/fatalerror.h"
  47 #include "gromacs/utility/futil.h"
  48 #include "gstat.h"
  49 #include "gromacs/legacyheaders/macros.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/units.h"
  52 #include "gromacs/utility/smalloc.h"
  53 #include "gromacs/commandline/pargs.h"
  54 #include "gromacs/legacyheaders/txtdump.h"
  55 #include "gromacs/legacyheaders/typedefs.h"
  56 #include "gromacs/math/vec.h"
  57 #include "gromacs/fileio/xvgr.h"
  58 #include "gromacs/legacyheaders/viewit.h"
  59 #include "correl.h"
  60 #include "gmx_ana.h"
  61 #include "gromacs/fft/fft.h"
  62 #include "gromacs/fileio/trxio.h"
  63
  64 enum {
  65     VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR
  66 };
  67
  68 static double FD(double Delta, double f)
  69 {
  70     return (2*pow(Delta, -4.5)*pow(f, 7.5) -
  71             6*pow(Delta, -3)*pow(f, 5) -
  72             pow(Delta, -1.5)*pow(f, 3.5) +
  73             6*pow(Delta, -1.5)*pow(f, 2.5) +
  74             2*f - 2);
  75 }
  76
  77 static double YYY(double f, double y)
  78 {
  79     return (2*pow(y*f, 3) - sqr(f)*y*(1+6*y) +
  80             (2+6*y)*f - 2);
  81 }
  82
  83 static double calc_compress(double y)
  84 {
  85     if (y == 1)
  86     {
  87         return 0;
  88     }
  89     return ((1+y+sqr(y)-pow(y, 3))/(pow(1-y, 3)));
  90 }
  91
  92 static double bisector(double Delta, double tol,
  93                        double ff0, double ff1,
  94                        double ff(double, double))
  95 {
  96     double fd0, fd, fd1, f, f0, f1;
  97     double tolmin = 1e-8;
  98
  99     f0 = ff0;
 100     f1 = ff1;
 101     if (tol < tolmin)
 102     {
 103         fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
 104         tol = tolmin;
 105     }
 106
 107     do
 108     {
 109         fd0 = ff(Delta, f0);
 110         fd1 = ff(Delta, f1);
 111         f   = (f0+f1)*0.5;
 112         fd  = ff(Delta, f);
 113         if (fd < 0)
 114         {
 115             f0 = f;
 116         }
 117         else if (fd > 0)
 118         {
 119             f1 = f;
 120         }
 121         else
 122         {
 123             return f;
 124         }
 125     }
 126     while ((f1-f0) > tol);
 127
 128     return f;
 129 }
 130
 131 static double calc_fluidicity(double Delta, double tol)
 132 {
 133     return bisector(Delta, tol, 0, 1, FD);
 134 }
 135
 136 static double calc_y(double f, double Delta, double toler)
 137 {
 138     double y1, y2;
 139
 140     y1 = pow(f/Delta, 1.5);
 141     y2 = bisector(f, toler, 0, 10000, YYY);
 142     if (fabs((y1-y2)/(y1+y2)) > 100*toler)
 143     {
 144         fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
 145                 y1, y2);
 146     }
 147
 148     return y1;
 149 }
 150
 151 static double calc_Shs(double f, double y)
 152 {
 153     double fy  = f*y;
 154
 155     return BOLTZ*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
 156 }
 157
 158 static real wCsolid(real nu, real beta)
 159 {
 160     real bhn = beta*PLANCK*nu;
 161     real ebn, koko;
 162
 163     if (bhn == 0)
 164     {
 165         return 1.0;
 166     }
 167     else
 168     {
 169         ebn  = exp(bhn);
 170         koko = sqr(1-ebn);
 171         return sqr(bhn)*ebn/koko;
 172     }
 173 }
 174
 175 static real wSsolid(real nu, real beta)
 176 {
 177     real bhn = beta*PLANCK*nu;
 178
 179     if (bhn == 0)
 180     {
 181         return 1;
 182     }
 183     else
 184     {
 185         return bhn/(exp(bhn)-1) - log(1-exp(-bhn));
 186     }
 187 }
 188
 189 static real wAsolid(real nu, real beta)
 190 {
 191     real bhn = beta*PLANCK*nu;
 192
 193     if (bhn == 0)
 194     {
 195         return 0;
 196     }
 197     else
 198     {
 199         return log((1-exp(-bhn))/(exp(-bhn/2))) - log(bhn);
 200     }
 201 }
 202
 203 static real wEsolid(real nu, real beta)
 204 {
 205     real bhn = beta*PLANCK*nu;
 206
 207     if (bhn == 0)
 208     {
 209         return 1;
 210     }
 211     else
 212     {
 213         return bhn/2 + bhn/(exp(bhn)-1)-1;
 214     }
 215 }
 216
 217 static void dump_fy(output_env_t oenv, real toler)
 218 {
 219     FILE       *fp;
 220     double      Delta, f, y, DD;
 221     const char *leg[] = { "f", "fy", "y" };
 222
 223     DD = pow(10.0, 0.125);
 224     fp = xvgropen("fy.xvg", "Fig. 2, Lin2003a", "Delta", "y or fy", oenv);
 225     xvgr_legend(fp, asize(leg), leg, oenv);
 226     if (output_env_get_print_xvgr_codes(oenv))
 227     {
 228         fprintf(fp, "@    world 1e-05, 0, 1000, 1\n");
 229         fprintf(fp, "@    xaxes scale Logarithmic\n");
 230     }
 231     for (Delta = 1e-5; (Delta <= 1000); Delta *= DD)
 232     {
 233         f = calc_fluidicity(Delta, toler);
 234         y = calc_y(f, Delta, toler);
 235         fprintf(fp, "%10g  %10g  %10g  %10g\n", Delta, f, f*y, y);
 236     }
 237     xvgrclose(fp);
 238 }
 239
 240 static void dump_w(output_env_t oenv, real beta)
 241 {
 242     FILE       *fp;
 243     double      nu;
 244     const char *leg[] = { "wCv", "wS", "wA", "wE" };
 245
 246     fp = xvgropen("w.xvg", "Fig. 1, Berens1983a", "\\f{12}b\\f{4}h\\f{12}n",
 247                   "w", oenv);
 248     xvgr_legend(fp, asize(leg), leg, oenv);
 249     for (nu = 1; (nu < 100); nu += 0.05)
 250     {
 251         fprintf(fp, "%10g  %10g  %10g  %10g  %10g\n", beta*PLANCK*nu,
 252                 wCsolid(nu, beta), wSsolid(nu, beta),
 253                 wAsolid(nu, beta), wEsolid(nu, beta));
 254     }
 255     xvgrclose(fp);
 256 }
 257
 258 int gmx_dos(int argc, char *argv[])
 259 {
 260     const char         *desc[] = {
 261         "[THISMODULE] computes the Density of States from a simulations.",
 262         "In order for this to be meaningful the velocities must be saved",
 263         "in the trajecotry with sufficiently high frequency such as to cover",
 264         "all vibrations. For flexible systems that would be around a few fs",
 265         "between saving. Properties based on the DoS are printed on the",
 266         "standard output."
 267     };
 268     const char         *bugs[] = {
 269         "This program needs a lot of memory: total usage equals the number of atoms times 3 times number of frames times 4 (or 8 when run in double precision)."
 270     };
 271     FILE               *fp, *fplog;
 272     t_topology          top;
 273     int                 ePBC = -1;
 274     t_trxframe          fr;
 275     matrix              box;
 276     int                 gnx;
 277     char                title[256];
 278     real                t0, t1, m;
 279     t_trxstatus        *status;
 280     int                 nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
 281     double              rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
 282     real              **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
 283     output_env_t        oenv;
 284     gmx_fft_t           fft;
 285     double              cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
 286     double              wCdiff, wSdiff, wAdiff, wEdiff;
 287
 288     static     gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalize = FALSE;
 289     static     gmx_bool bRecip   = FALSE, bDump = FALSE;
 290     static     real     Temp     = 298.15, toler = 1e-6;
 291     t_pargs             pa[]     = {
 292         { "-v", FALSE, etBOOL, {&bVerbose},
 293           "Be loud and noisy." },
 294         { "-recip", FALSE, etBOOL, {&bRecip},
 295           "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
 296         { "-abs", FALSE, etBOOL, {&bAbsolute},
 297           "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
 298         { "-normdos", FALSE, etBOOL, {&bNormalize},
 299           "Normalize the DoS such that it adds up to 3N. This is a hack that should not be necessary." },
 300         { "-T", FALSE, etREAL, {&Temp},
 301           "Temperature in the simulation" },
 302         { "-toler", FALSE, etREAL, {&toler},
 303           "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" },
 304         { "-dump", FALSE, etBOOL, {&bDump},
 305           "[HIDDEN]Dump the y/fy plot corresponding to Fig. 2 inLin2003a and the and the weighting functions corresponding to Fig. 1 in Berens1983a." }
 306     };
 307
 308     t_filenm            fnm[] = {
 309         { efTRN, "-f",    NULL,    ffREAD  },
 310         { efTPX, "-s",    NULL,    ffREAD  },
 311         { efNDX, NULL,    NULL,    ffOPTRD },
 312         { efXVG, "-vacf", "vacf",  ffWRITE },
 313         { efXVG, "-mvacf", "mvacf", ffWRITE },
 314         { efXVG, "-dos",  "dos",   ffWRITE },
 315         { efLOG, "-g",    "dos",   ffWRITE },
 316     };
 317 #define NFILE asize(fnm)
 318     int                 npargs;
 319     t_pargs            *ppa;
 320     const char         *DoSlegend[] = {
 321         "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
 322     };
 323
 324     npargs = asize(pa);
 325     ppa    = add_acf_pargs(&npargs, pa);
 326     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 327                            NFILE, fnm, npargs, ppa, asize(desc), desc,
 328                            asize(bugs), bugs, &oenv))
 329     {
 330         return 0;
 331     }
 332
 333     beta = 1/(Temp*BOLTZ);
 334     if (bDump)
 335     {
 336         printf("Dumping reference figures. Thanks for your patience.\n");
 337         dump_fy(oenv, toler);
 338         dump_w(oenv, beta);
 339         exit(0);
 340     }
 341
 342     fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
 343     fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
 344     please_cite(fplog, "Pascal2011a");
 345     please_cite(fplog, "Caleman2011b");
 346
 347     read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
 348                   TRUE);
 349     V     = det(box);
 350     tmass = 0;
 351     for (i = 0; (i < top.atoms.nr); i++)
 352     {
 353         tmass += top.atoms.atom[i].m;
 354     }
 355
 356     Natom = top.atoms.nr;
 357     Nmol  = top.mols.nr;
 358     gnx   = Natom*DIM;
 359
 360     /* Correlation stuff */
 361     snew(c1, gnx);
 362     for (i = 0; (i < gnx); i++)
 363     {
 364         c1[i] = NULL;
 365     }
 366
 367     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
 368     t0 = fr.time;
 369
 370     n_alloc = 0;
 371     nframes = 0;
 372     Vsum    = V2sum = 0;
 373     nV      = 0;
 374     do
 375     {
 376         if (fr.bBox)
 377         {
 378             V      = det(fr.box);
 379             V2sum += V*V;
 380             Vsum  += V;
 381             nV++;
 382         }
 383         if (nframes >= n_alloc)
 384         {
 385             n_alloc += 100;
 386             for (i = 0; i < gnx; i++)
 387             {
 388                 srenew(c1[i], n_alloc);
 389             }
 390         }
 391         for (i = 0; i < gnx; i += DIM)
 392         {
 393             c1[i+XX][nframes] = fr.v[i/DIM][XX];
 394             c1[i+YY][nframes] = fr.v[i/DIM][YY];
 395             c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
 396         }
 397
 398         t1 = fr.time;
 399
 400         nframes++;
 401     }
 402     while (read_next_frame(oenv, status, &fr));
 403
 404     close_trj(status);
 405
 406     dt = (t1-t0)/(nframes-1);
 407     if (nV > 0)
 408     {
 409         V = Vsum/nV;
 410     }
 411     if (bVerbose)
 412     {
 413         printf("Going to do %d fourier transforms of length %d. Hang on.\n",
 414                gnx, nframes);
 415     }
 416     low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
 417                     FALSE, FALSE, -1, -1, 0);
 418     snew(dos, DOS_NR);
 419     for (j = 0; (j < DOS_NR); j++)
 420     {
 421         snew(dos[j], nframes+4);
 422     }
 423
 424     if (bVerbose)
 425     {
 426         printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
 427     }
 428     for (i = 0; (i < gnx); i += DIM)
 429     {
 430         mi = top.atoms.atom[i/DIM].m;
 431         for (j = 0; (j < nframes/2); j++)
 432         {
 433             c1j            = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
 434             dos[VACF][j]  += c1j/Natom;
 435             dos[MVACF][j] += mi*c1j;
 436         }
 437     }
 438     fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF",
 439                   "Time (ps)", "C(t)", oenv);
 440     snew(tt, nframes/2);
 441     for (j = 0; (j < nframes/2); j++)
 442     {
 443         tt[j] = j*dt;
 444         fprintf(fp, "%10g  %10g\n", tt[j], dos[VACF][j]);
 445     }
 446     xvgrclose(fp);
 447     fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF",
 448                   "Time (ps)", "C(t)", oenv);
 449     for (j = 0; (j < nframes/2); j++)
 450     {
 451         fprintf(fp, "%10g  %10g\n", tt[j], dos[MVACF][j]);
 452     }
 453     xvgrclose(fp);
 454
 455     if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
 456                                         GMX_FFT_FLAG_NONE)) != 0)
 457     {
 458         gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
 459     }
 460     if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
 461                                    (void *)dos[MVACF], (void *)dos[DOS])) != 0)
 462     {
 463         gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
 464     }
 465
 466     /* First compute the DoS */
 467     /* Magic factor of 8 included now. */
 468     bfac = 8*dt*beta/2;
 469     dos2 = 0;
 470     snew(nu, nframes/4);
 471     for (j = 0; (j < nframes/4); j++)
 472     {
 473         nu[j] = 2*j/(t1-t0);
 474         dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
 475         if (bAbsolute)
 476         {
 477             dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
 478         }
 479         else
 480         {
 481             dos[DOS][j] = bfac*dos[DOS][2*j];
 482         }
 483     }
 484     /* Normalize it */
 485     dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
 486     if (bNormalize)
 487     {
 488         for (j = 0; (j < nframes/4); j++)
 489         {
 490             dos[DOS][j] *= 3*Natom/dostot;
 491         }
 492     }
 493
 494     /* Now analyze it */
 495     DoS0 = dos[DOS][0];
 496
 497     /* Note this eqn. is incorrect in Pascal2011a! */
 498     Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
 499              pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
 500     f     = calc_fluidicity(Delta, toler);
 501     y     = calc_y(f, Delta, toler);
 502     z     = calc_compress(y);
 503     Sig   = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
 504     Shs   = Sig+calc_Shs(f, y);
 505     rho   = (tmass*AMU)/(V*NANO*NANO*NANO);
 506     sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
 507
 508     fprintf(fplog, "System = \"%s\"\n", title);
 509     fprintf(fplog, "Nmol = %d\n", Nmol);
 510     fprintf(fplog, "Natom = %d\n", Natom);
 511     fprintf(fplog, "dt = %g ps\n", dt);
 512     fprintf(fplog, "tmass = %g amu\n", tmass);
 513     fprintf(fplog, "V = %g nm^3\n", V);
 514     fprintf(fplog, "rho = %g g/l\n", rho);
 515     fprintf(fplog, "T = %g K\n", Temp);
 516     fprintf(fplog, "beta = %g mol/kJ\n", beta);
 517
 518     fprintf(fplog, "\nDoS parameters\n");
 519     fprintf(fplog, "Delta = %g\n", Delta);
 520     fprintf(fplog, "fluidicity = %g\n", f);
 521     fprintf(fplog, "hard sphere packing fraction = %g\n", y);
 522     fprintf(fplog, "hard sphere compressibility = %g\n", z);
 523     fprintf(fplog, "ideal gas entropy = %g\n", Sig);
 524     fprintf(fplog, "hard sphere entropy = %g\n", Shs);
 525     fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
 526     fprintf(fplog, "DoS0 = %g\n", DoS0);
 527     fprintf(fplog, "Dos2 = %g\n", dos2);
 528     fprintf(fplog, "DoSTot = %g\n", dostot);
 529
 530     /* Now compute solid (2) and diffusive (3) components */
 531     fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
 532                   bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
 533                   "\\f{4}S(\\f{12}n\\f{4})", oenv);
 534     xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
 535     recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
 536     for (j = 0; (j < nframes/4); j++)
 537     {
 538         dos[DOS_DIFF][j]  = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
 539         dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
 540         fprintf(fp, "%10g  %10g  %10g  %10g\n",
 541                 recip_fac*nu[j],
 542                 dos[DOS][j]/recip_fac,
 543                 dos[DOS_SOLID][j]/recip_fac,
 544                 dos[DOS_DIFF][j]/recip_fac);
 545     }
 546     xvgrclose(fp);
 547
 548     /* Finally analyze the results! */
 549     wCdiff = 0.5;
 550     wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
 551     wEdiff = 0.5;
 552     wAdiff = wEdiff-wSdiff;
 553     for (j = 0; (j < nframes/4); j++)
 554     {
 555         dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
 556                           dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
 557         dos[DOS_S][j]  = (dos[DOS_DIFF][j]*wSdiff +
 558                           dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
 559         dos[DOS_A][j]  = (dos[DOS_DIFF][j]*wAdiff +
 560                           dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
 561         dos[DOS_E][j]  = (dos[DOS_DIFF][j]*wEdiff +
 562                           dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
 563     }
 564     DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
 565     DiffCoeff = 1000*DiffCoeff/3.0;
 566     fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
 567             DiffCoeff);
 568     fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
 569             1000*DoS0/(12*tmass*beta));
 570
 571     cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
 572                                    nframes/4, &stddev);
 573     fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
 574
 575     /*
 576        S  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
 577                                    nframes/4,&stddev);
 578        fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
 579        A  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,
 580                                    nframes/4,&stddev);
 581        fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol);
 582        E  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,
 583                                    nframes/4,&stddev);
 584        fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol);
 585      */
 586     fprintf(fplog, "\nArrivederci!\n");
 587     gmx_fio_fclose(fplog);
 588
 589     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");
 590
 591     return 0;
 592 }