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