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