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