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 #include <stdio.h>
  39 #include <math.h>
  40
  41 #include "gromacs/fileio/confio.h"
  42 #include "copyrite.h"
  43 #include "gmx_fatal.h"
  44 #include "gromacs/fileio/futil.h"
  45 #include "gstat.h"
  46 #include "macros.h"
  47 #include "gromacs/math/utilities.h"
  48 #include "physics.h"
  49 #include "index.h"
  50 #include "smalloc.h"
  51 #include "gromacs/commandline/pargs.h"
  52 #include <string.h>
  53 #include "sysstuff.h"
  54 #include "txtdump.h"
  55 #include "typedefs.h"
  56 #include "vec.h"
  57 #include "strdb.h"
  58 #include "xvgr.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     fprintf(fp, "@    world 1e-05, 0, 1000, 1\n");
 227     fprintf(fp, "@    xaxes scale Logarithmic\n");
 228     for (Delta = 1e-5; (Delta <= 1000); Delta *= DD)
 229     {
 230         f = calc_fluidicity(Delta, toler);
 231         y = calc_y(f, Delta, toler);
 232         fprintf(fp, "%10g  %10g  %10g  %10g\n", Delta, f, f*y, y);
 233     }
 234     xvgrclose(fp);
 235 }
 236
 237 static void dump_w(output_env_t oenv, real beta)
 238 {
 239     FILE       *fp;
 240     double      nu;
 241     const char *leg[] = { "wCv", "wS", "wA", "wE" };
 242
 243     fp = xvgropen("w.xvg", "Fig. 1, Berens1983a", "\\f{12}b\\f{4}h\\f{12}n",
 244                   "w", oenv);
 245     xvgr_legend(fp, asize(leg), leg, oenv);
 246     for (nu = 1; (nu < 100); nu += 0.05)
 247     {
 248         fprintf(fp, "%10g  %10g  %10g  %10g  %10g\n", beta*PLANCK*nu,
 249                 wCsolid(nu, beta), wSsolid(nu, beta),
 250                 wAsolid(nu, beta), wEsolid(nu, beta));
 251     }
 252     xvgrclose(fp);
 253 }
 254
 255 int gmx_dos(int argc, char *argv[])
 256 {
 257     const char         *desc[] = {
 258         "[THISMODULE] computes the Density of States from a simulations.",
 259         "In order for this to be meaningful the velocities must be saved",
 260         "in the trajecotry with sufficiently high frequency such as to cover",
 261         "all vibrations. For flexible systems that would be around a few fs",
 262         "between saving. Properties based on the DoS are printed on the",
 263         "standard output."
 264     };
 265     const char         *bugs[] = {
 266         "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)."
 267     };
 268     FILE               *fp, *fplog;
 269     t_topology          top;
 270     int                 ePBC = -1;
 271     t_trxframe          fr;
 272     matrix              box;
 273     int                 gnx;
 274     char                title[256];
 275     real                t0, t1, m;
 276     t_trxstatus        *status;
 277     int                 nV, nframes, n_alloc, i, j, k, l, fftcode, Nmol, Natom;
 278     double              rho, dt, V2sum, Vsum, V, tmass, dostot, dos2, dosabs;
 279     real              **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
 280     output_env_t        oenv;
 281     gmx_fft_t           fft;
 282     double              cP, S, A, E, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
 283     double              wCdiff, wSdiff, wAdiff, wEdiff;
 284
 285     static     gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalize = FALSE;
 286     static     gmx_bool bRecip   = FALSE, bDump = FALSE;
 287     static     real     Temp     = 298.15, toler = 1e-6;
 288     t_pargs             pa[]     = {
 289         { "-v", FALSE, etBOOL, {&bVerbose},
 290           "Be loud and noisy." },
 291         { "-recip", FALSE, etBOOL, {&bRecip},
 292           "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
 293         { "-abs", FALSE, etBOOL, {&bAbsolute},
 294           "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
 295         { "-normdos", FALSE, etBOOL, {&bNormalize},
 296           "Normalize the DoS such that it adds up to 3N. This is a hack that should not be necessary." },
 297         { "-T", FALSE, etREAL, {&Temp},
 298           "Temperature in the simulation" },
 299         { "-toler", FALSE, etREAL, {&toler},
 300           "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" },
 301         { "-dump", FALSE, etBOOL, {&bDump},
 302           "[HIDDEN]Dump the y/fy plot corresponding to Fig. 2 inLin2003a and the and the weighting functions corresponding to Fig. 1 in Berens1983a." }
 303     };
 304
 305     t_filenm            fnm[] = {
 306         { efTRN, "-f",    NULL,    ffREAD  },
 307         { efTPX, "-s",    NULL,    ffREAD  },
 308         { efNDX, NULL,    NULL,    ffOPTRD },
 309         { efXVG, "-vacf", "vacf",  ffWRITE },
 310         { efXVG, "-mvacf", "mvacf", ffWRITE },
 311         { efXVG, "-dos",  "dos",   ffWRITE },
 312         { efLOG, "-g",    "dos",   ffWRITE },
 313     };
 314 #define NFILE asize(fnm)
 315     int                 npargs;
 316     t_pargs            *ppa;
 317     const char         *DoSlegend[] = {
 318         "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
 319     };
 320
 321     npargs = asize(pa);
 322     ppa    = add_acf_pargs(&npargs, pa);
 323     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 324                            NFILE, fnm, npargs, ppa, asize(desc), desc,
 325                            asize(bugs), bugs, &oenv))
 326     {
 327         return 0;
 328     }
 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);
 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     return 0;
 589 }