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     if (output_env_get_print_xvgr_codes(oenv))
 226     {
 227         fprintf(fp, "@    world 1e-05, 0, 1000, 1\n");
 228         fprintf(fp, "@    xaxes scale Logarithmic\n");
 229     }
 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         "[THISMODULE] 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     npargs = asize(pa);
 324     ppa    = add_acf_pargs(&npargs, pa);
 325     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 326                            NFILE, fnm, npargs, ppa, asize(desc), desc,
 327                            asize(bugs), bugs, &oenv))
 328     {
 329         return 0;
 330     }
 331
 332     beta = 1/(Temp*BOLTZ);
 333     if (bDump)
 334     {
 335         printf("Dumping reference figures. Thanks for your patience.\n");
 336         dump_fy(oenv, toler);
 337         dump_w(oenv, beta);
 338         exit(0);
 339     }
 340
 341     fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
 342     fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
 343     please_cite(fplog, "Pascal2011a");
 344     please_cite(fplog, "Caleman2011b");
 345
 346     read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
 347                   TRUE);
 348     V     = det(box);
 349     tmass = 0;
 350     for (i = 0; (i < top.atoms.nr); i++)
 351     {
 352         tmass += top.atoms.atom[i].m;
 353     }
 354
 355     Natom = top.atoms.nr;
 356     Nmol  = top.mols.nr;
 357     gnx   = Natom*DIM;
 358
 359     /* Correlation stuff */
 360     snew(c1, gnx);
 361     for (i = 0; (i < gnx); i++)
 362     {
 363         c1[i] = NULL;
 364     }
 365
 366     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
 367     t0 = fr.time;
 368
 369     n_alloc = 0;
 370     nframes = 0;
 371     Vsum    = V2sum = 0;
 372     nV      = 0;
 373     do
 374     {
 375         if (fr.bBox)
 376         {
 377             V      = det(fr.box);
 378             V2sum += V*V;
 379             Vsum  += V;
 380             nV++;
 381         }
 382         if (nframes >= n_alloc)
 383         {
 384             n_alloc += 100;
 385             for (i = 0; i < gnx; i++)
 386             {
 387                 srenew(c1[i], n_alloc);
 388             }
 389         }
 390         for (i = 0; i < gnx; i += DIM)
 391         {
 392             c1[i+XX][nframes] = fr.v[i/DIM][XX];
 393             c1[i+YY][nframes] = fr.v[i/DIM][YY];
 394             c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
 395         }
 396
 397         t1 = fr.time;
 398
 399         nframes++;
 400     }
 401     while (read_next_frame(oenv, status, &fr));
 402
 403     close_trj(status);
 404
 405     dt = (t1-t0)/(nframes-1);
 406     if (nV > 0)
 407     {
 408         V = Vsum/nV;
 409     }
 410     if (bVerbose)
 411     {
 412         printf("Going to do %d fourier transforms of length %d. Hang on.\n",
 413                gnx, nframes);
 414     }
 415     low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
 416                     FALSE, FALSE, -1, -1, 0);
 417     snew(dos, DOS_NR);
 418     for (j = 0; (j < DOS_NR); j++)
 419     {
 420         snew(dos[j], nframes+4);
 421     }
 422
 423     if (bVerbose)
 424     {
 425         printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
 426     }
 427     for (i = 0; (i < gnx); i += DIM)
 428     {
 429         mi = top.atoms.atom[i/DIM].m;
 430         for (j = 0; (j < nframes/2); j++)
 431         {
 432             c1j            = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
 433             dos[VACF][j]  += c1j/Natom;
 434             dos[MVACF][j] += mi*c1j;
 435         }
 436     }
 437     fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF",
 438                   "Time (ps)", "C(t)", oenv);
 439     snew(tt, nframes/2);
 440     for (j = 0; (j < nframes/2); j++)
 441     {
 442         tt[j] = j*dt;
 443         fprintf(fp, "%10g  %10g\n", tt[j], dos[VACF][j]);
 444     }
 445     xvgrclose(fp);
 446     fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF",
 447                   "Time (ps)", "C(t)", oenv);
 448     for (j = 0; (j < nframes/2); j++)
 449     {
 450         fprintf(fp, "%10g  %10g\n", tt[j], dos[MVACF][j]);
 451     }
 452     xvgrclose(fp);
 453
 454     if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
 455                                         GMX_FFT_FLAG_NONE)) != 0)
 456     {
 457         gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
 458     }
 459     if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
 460                                    (void *)dos[MVACF], (void *)dos[DOS])) != 0)
 461     {
 462         gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
 463     }
 464
 465     /* First compute the DoS */
 466     /* Magic factor of 8 included now. */
 467     bfac = 8*dt*beta/2;
 468     dos2 = 0;
 469     snew(nu, nframes/4);
 470     for (j = 0; (j < nframes/4); j++)
 471     {
 472         nu[j] = 2*j/(t1-t0);
 473         dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
 474         if (bAbsolute)
 475         {
 476             dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
 477         }
 478         else
 479         {
 480             dos[DOS][j] = bfac*dos[DOS][2*j];
 481         }
 482     }
 483     /* Normalize it */
 484     dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
 485     if (bNormalize)
 486     {
 487         for (j = 0; (j < nframes/4); j++)
 488         {
 489             dos[DOS][j] *= 3*Natom/dostot;
 490         }
 491     }
 492
 493     /* Now analyze it */
 494     DoS0 = dos[DOS][0];
 495
 496     /* Note this eqn. is incorrect in Pascal2011a! */
 497     Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
 498              pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
 499     f     = calc_fluidicity(Delta, toler);
 500     y     = calc_y(f, Delta, toler);
 501     z     = calc_compress(y);
 502     Sig   = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
 503     Shs   = Sig+calc_Shs(f, y);
 504     rho   = (tmass*AMU)/(V*NANO*NANO*NANO);
 505     sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
 506
 507     fprintf(fplog, "System = \"%s\"\n", title);
 508     fprintf(fplog, "Nmol = %d\n", Nmol);
 509     fprintf(fplog, "Natom = %d\n", Natom);
 510     fprintf(fplog, "dt = %g ps\n", dt);
 511     fprintf(fplog, "tmass = %g amu\n", tmass);
 512     fprintf(fplog, "V = %g nm^3\n", V);
 513     fprintf(fplog, "rho = %g g/l\n", rho);
 514     fprintf(fplog, "T = %g K\n", Temp);
 515     fprintf(fplog, "beta = %g mol/kJ\n", beta);
 516
 517     fprintf(fplog, "\nDoS parameters\n");
 518     fprintf(fplog, "Delta = %g\n", Delta);
 519     fprintf(fplog, "fluidicity = %g\n", f);
 520     fprintf(fplog, "hard sphere packing fraction = %g\n", y);
 521     fprintf(fplog, "hard sphere compressibility = %g\n", z);
 522     fprintf(fplog, "ideal gas entropy = %g\n", Sig);
 523     fprintf(fplog, "hard sphere entropy = %g\n", Shs);
 524     fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
 525     fprintf(fplog, "DoS0 = %g\n", DoS0);
 526     fprintf(fplog, "Dos2 = %g\n", dos2);
 527     fprintf(fplog, "DoSTot = %g\n", dostot);
 528
 529     /* Now compute solid (2) and diffusive (3) components */
 530     fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
 531                   bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
 532                   "\\f{4}S(\\f{12}n\\f{4})", oenv);
 533     xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
 534     recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
 535     for (j = 0; (j < nframes/4); j++)
 536     {
 537         dos[DOS_DIFF][j]  = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
 538         dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
 539         fprintf(fp, "%10g  %10g  %10g  %10g\n",
 540                 recip_fac*nu[j],
 541                 dos[DOS][j]/recip_fac,
 542                 dos[DOS_SOLID][j]/recip_fac,
 543                 dos[DOS_DIFF][j]/recip_fac);
 544     }
 545     xvgrclose(fp);
 546
 547     /* Finally analyze the results! */
 548     wCdiff = 0.5;
 549     wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
 550     wEdiff = 0.5;
 551     wAdiff = wEdiff-wSdiff;
 552     for (j = 0; (j < nframes/4); j++)
 553     {
 554         dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
 555                           dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
 556         dos[DOS_S][j]  = (dos[DOS_DIFF][j]*wSdiff +
 557                           dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
 558         dos[DOS_A][j]  = (dos[DOS_DIFF][j]*wAdiff +
 559                           dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
 560         dos[DOS_E][j]  = (dos[DOS_DIFF][j]*wEdiff +
 561                           dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
 562     }
 563     DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
 564     DiffCoeff = 1000*DiffCoeff/3.0;
 565     fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
 566             DiffCoeff);
 567     fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
 568             1000*DoS0/(12*tmass*beta));
 569
 570     cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
 571                                    nframes/4, &stddev);
 572     fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
 573
 574     /*
 575        S  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
 576                                    nframes/4,&stddev);
 577        fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
 578        A  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,
 579                                    nframes/4,&stddev);
 580        fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol);
 581        E  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,
 582                                    nframes/4,&stddev);
 583        fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol);
 584      */
 585     fprintf(fplog, "\nArrivederci!\n");
 586     gmx_fio_fclose(fplog);
 587
 588     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");
 589
 590     return 0;
 591 }