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 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  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 "copyrite.h"
  46 #include "gromacs/utility/fatalerror.h"
  47 #include "gromacs/utility/futil.h"
  48 #include "gstat.h"
  49 #include "macros.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/units.h"
  52 #include "index.h"
  53 #include "gromacs/utility/smalloc.h"
  54 #include "gromacs/commandline/pargs.h"
  55 #include "txtdump.h"
  56 #include "typedefs.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/fileio/xvgr.h"
  59 #include "viewit.h"
  60 #include "correl.h"
  61 #include "gmx_ana.h"
  62 #include "gromacs/fft/fft.h"
  63 #include "gromacs/fileio/trxio.h"
  64
  65 enum {
  66     VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR
  67 };
  68
  69 static double FD(double Delta, double f)
  70 {
  71     return (2*pow(Delta, -4.5)*pow(f, 7.5) -
  72             6*pow(Delta, -3)*pow(f, 5) -
  73             pow(Delta, -1.5)*pow(f, 3.5) +
  74             6*pow(Delta, -1.5)*pow(f, 2.5) +
  75             2*f - 2);
  76 }
  77
  78 static double YYY(double f, double y)
  79 {
  80     return (2*pow(y*f, 3) - sqr(f)*y*(1+6*y) +
  81             (2+6*y)*f - 2);
  82 }
  83
  84 static double calc_compress(double y)
  85 {
  86     if (y == 1)
  87     {
  88         return 0;
  89     }
  90     return ((1+y+sqr(y)-pow(y, 3))/(pow(1-y, 3)));
  91 }
  92
  93 static double bisector(double Delta, double tol,
  94                        double ff0, double ff1,
  95                        double ff(double, double))
  96 {
  97     double fd0, fd, fd1, f, f0, f1;
  98     double tolmin = 1e-8;
  99
 100     f0 = ff0;
 101     f1 = ff1;
 102     if (tol < tolmin)
 103     {
 104         fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
 105         tol = tolmin;
 106     }
 107
 108     do
 109     {
 110         fd0 = ff(Delta, f0);
 111         fd1 = ff(Delta, f1);
 112         f   = (f0+f1)*0.5;
 113         fd  = ff(Delta, f);
 114         if (fd < 0)
 115         {
 116             f0 = f;
 117         }
 118         else if (fd > 0)
 119         {
 120             f1 = f;
 121         }
 122         else
 123         {
 124             return f;
 125         }
 126     }
 127     while ((f1-f0) > tol);
 128
 129     return f;
 130 }
 131
 132 static double calc_fluidicity(double Delta, double tol)
 133 {
 134     return bisector(Delta, tol, 0, 1, FD);
 135 }
 136
 137 static double calc_y(double f, double Delta, double toler)
 138 {
 139     double y1, y2;
 140
 141     y1 = pow(f/Delta, 1.5);
 142     y2 = bisector(f, toler, 0, 10000, YYY);
 143     if (fabs((y1-y2)/(y1+y2)) > 100*toler)
 144     {
 145         fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
 146                 y1, y2);
 147     }
 148
 149     return y1;
 150 }
 151
 152 static double calc_Shs(double f, double y)
 153 {
 154     double fy  = f*y;
 155
 156     return BOLTZ*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
 157 }
 158
 159 static real wCsolid(real nu, real beta)
 160 {
 161     real bhn = beta*PLANCK*nu;
 162     real ebn, koko;
 163
 164     if (bhn == 0)
 165     {
 166         return 1.0;
 167     }
 168     else
 169     {
 170         ebn  = exp(bhn);
 171         koko = sqr(1-ebn);
 172         return sqr(bhn)*ebn/koko;
 173     }
 174 }
 175
 176 static real wSsolid(real nu, real beta)
 177 {
 178     real bhn = beta*PLANCK*nu;
 179
 180     if (bhn == 0)
 181     {
 182         return 1;
 183     }
 184     else
 185     {
 186         return bhn/(exp(bhn)-1) - log(1-exp(-bhn));
 187     }
 188 }
 189
 190 static real wAsolid(real nu, real beta)
 191 {
 192     real bhn = beta*PLANCK*nu;
 193
 194     if (bhn == 0)
 195     {
 196         return 0;
 197     }
 198     else
 199     {
 200         return log((1-exp(-bhn))/(exp(-bhn/2))) - log(bhn);
 201     }
 202 }
 203
 204 static real wEsolid(real nu, real beta)
 205 {
 206     real bhn = beta*PLANCK*nu;
 207
 208     if (bhn == 0)
 209     {
 210         return 1;
 211     }
 212     else
 213     {
 214         return bhn/2 + bhn/(exp(bhn)-1)-1;
 215     }
 216 }
 217
 218 static void dump_fy(output_env_t oenv, real toler)
 219 {
 220     FILE       *fp;
 221     double      Delta, f, y, DD;
 222     const char *leg[] = { "f", "fy", "y" };
 223
 224     DD = pow(10.0, 0.125);
 225     fp = xvgropen("fy.xvg", "Fig. 2, Lin2003a", "Delta", "y or fy", oenv);
 226     xvgr_legend(fp, asize(leg), leg, oenv);
 227     fprintf(fp, "@    world 1e-05, 0, 1000, 1\n");
 228     fprintf(fp, "@    xaxes scale Logarithmic\n");
 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 | PCA_BE_NICE,
 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 }