src/gromacs/gmxana/gmx_dos.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <stdio.h>
  39 #include <math.h>
  40
  41 #include "confio.h"
  42 #include "copyrite.h"
  43 #include "gmx_fatal.h"
  44 #include "futil.h"
  45 #include "gstat.h"
  46 #include "macros.h"
  47 #include "maths.h"
  48 #include "physics.h"
  49 #include "index.h"
  50 #include "smalloc.h"
  51 #include "statutil.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
  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         "[TT]g_dos[tt] 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     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     beta = 1/(Temp*BOLTZ);
 327     if (bDump)
 328     {
 329         printf("Dumping reference figures. Thanks for your patience.\n");
 330         dump_fy(oenv, toler);
 331         dump_w(oenv, beta);
 332         exit(0);
 333     }
 334
 335     fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
 336     fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
 337     please_cite(fplog, "Pascal2011a");
 338     please_cite(fplog, "Caleman2011b");
 339
 340     read_tps_conf(ftp2fn(efTPX, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box,
 341                   TRUE);
 342     V     = det(box);
 343     tmass = 0;
 344     for (i = 0; (i < top.atoms.nr); i++)
 345     {
 346         tmass += top.atoms.atom[i].m;
 347     }
 348
 349     Natom = top.atoms.nr;
 350     Nmol  = top.mols.nr;
 351     gnx   = Natom*DIM;
 352
 353     /* Correlation stuff */
 354     snew(c1, gnx);
 355     for (i = 0; (i < gnx); i++)
 356     {
 357         c1[i] = NULL;
 358     }
 359
 360     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
 361     t0 = fr.time;
 362
 363     n_alloc = 0;
 364     nframes = 0;
 365     Vsum    = V2sum = 0;
 366     nV      = 0;
 367     do
 368     {
 369         if (fr.bBox)
 370         {
 371             V      = det(fr.box);
 372             V2sum += V*V;
 373             Vsum  += V;
 374             nV++;
 375         }
 376         if (nframes >= n_alloc)
 377         {
 378             n_alloc += 100;
 379             for (i = 0; i < gnx; i++)
 380             {
 381                 srenew(c1[i], n_alloc);
 382             }
 383         }
 384         for (i = 0; i < gnx; i += DIM)
 385         {
 386             c1[i+XX][nframes] = fr.v[i/DIM][XX];
 387             c1[i+YY][nframes] = fr.v[i/DIM][YY];
 388             c1[i+ZZ][nframes] = fr.v[i/DIM][ZZ];
 389         }
 390
 391         t1 = fr.time;
 392
 393         nframes++;
 394     }
 395     while (read_next_frame(oenv, status, &fr));
 396
 397     close_trj(status);
 398
 399     dt = (t1-t0)/(nframes-1);
 400     if (nV > 0)
 401     {
 402         V = Vsum/nV;
 403     }
 404     if (bVerbose)
 405     {
 406         printf("Going to do %d fourier transforms of length %d. Hang on.\n",
 407                gnx, nframes);
 408     }
 409     low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
 410                     FALSE, FALSE, -1, -1, 0, 0);
 411     snew(dos, DOS_NR);
 412     for (j = 0; (j < DOS_NR); j++)
 413     {
 414         snew(dos[j], nframes+4);
 415     }
 416
 417     if (bVerbose)
 418     {
 419         printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
 420     }
 421     for (i = 0; (i < gnx); i += DIM)
 422     {
 423         mi = top.atoms.atom[i/DIM].m;
 424         for (j = 0; (j < nframes/2); j++)
 425         {
 426             c1j            = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
 427             dos[VACF][j]  += c1j/Natom;
 428             dos[MVACF][j] += mi*c1j;
 429         }
 430     }
 431     fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity ACF",
 432                   "Time (ps)", "C(t)", oenv);
 433     snew(tt, nframes/2);
 434     for (j = 0; (j < nframes/2); j++)
 435     {
 436         tt[j] = j*dt;
 437         fprintf(fp, "%10g  %10g\n", tt[j], dos[VACF][j]);
 438     }
 439     xvgrclose(fp);
 440     fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity ACF",
 441                   "Time (ps)", "C(t)", oenv);
 442     for (j = 0; (j < nframes/2); j++)
 443     {
 444         fprintf(fp, "%10g  %10g\n", tt[j], dos[MVACF][j]);
 445     }
 446     xvgrclose(fp);
 447
 448     if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
 449                                         GMX_FFT_FLAG_NONE)) != 0)
 450     {
 451         gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
 452     }
 453     if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
 454                                    (void *)dos[MVACF], (void *)dos[DOS])) != 0)
 455     {
 456         gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
 457     }
 458
 459     /* First compute the DoS */
 460     /* Magic factor of 8 included now. */
 461     bfac = 8*dt*beta/2;
 462     dos2 = 0;
 463     snew(nu, nframes/4);
 464     for (j = 0; (j < nframes/4); j++)
 465     {
 466         nu[j] = 2*j/(t1-t0);
 467         dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
 468         if (bAbsolute)
 469         {
 470             dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
 471         }
 472         else
 473         {
 474             dos[DOS][j] = bfac*dos[DOS][2*j];
 475         }
 476     }
 477     /* Normalize it */
 478     dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
 479     if (bNormalize)
 480     {
 481         for (j = 0; (j < nframes/4); j++)
 482         {
 483             dos[DOS][j] *= 3*Natom/dostot;
 484         }
 485     }
 486
 487     /* Now analyze it */
 488     DoS0 = dos[DOS][0];
 489
 490     /* Note this eqn. is incorrect in Pascal2011a! */
 491     Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
 492              pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
 493     f     = calc_fluidicity(Delta, toler);
 494     y     = calc_y(f, Delta, toler);
 495     z     = calc_compress(y);
 496     Sig   = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
 497     Shs   = Sig+calc_Shs(f, y);
 498     rho   = (tmass*AMU)/(V*NANO*NANO*NANO);
 499     sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
 500
 501     fprintf(fplog, "System = \"%s\"\n", title);
 502     fprintf(fplog, "Nmol = %d\n", Nmol);
 503     fprintf(fplog, "Natom = %d\n", Natom);
 504     fprintf(fplog, "dt = %g ps\n", dt);
 505     fprintf(fplog, "tmass = %g amu\n", tmass);
 506     fprintf(fplog, "V = %g nm^3\n", V);
 507     fprintf(fplog, "rho = %g g/l\n", rho);
 508     fprintf(fplog, "T = %g K\n", Temp);
 509     fprintf(fplog, "beta = %g mol/kJ\n", beta);
 510
 511     fprintf(fplog, "\nDoS parameters\n");
 512     fprintf(fplog, "Delta = %g\n", Delta);
 513     fprintf(fplog, "fluidicity = %g\n", f);
 514     fprintf(fplog, "hard sphere packing fraction = %g\n", y);
 515     fprintf(fplog, "hard sphere compressibility = %g\n", z);
 516     fprintf(fplog, "ideal gas entropy = %g\n", Sig);
 517     fprintf(fplog, "hard sphere entropy = %g\n", Shs);
 518     fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
 519     fprintf(fplog, "DoS0 = %g\n", DoS0);
 520     fprintf(fplog, "Dos2 = %g\n", dos2);
 521     fprintf(fplog, "DoSTot = %g\n", dostot);
 522
 523     /* Now compute solid (2) and diffusive (3) components */
 524     fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
 525                   bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
 526                   "\\f{4}S(\\f{12}n\\f{4})", oenv);
 527     xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
 528     recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
 529     for (j = 0; (j < nframes/4); j++)
 530     {
 531         dos[DOS_DIFF][j]  = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
 532         dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
 533         fprintf(fp, "%10g  %10g  %10g  %10g\n",
 534                 recip_fac*nu[j],
 535                 dos[DOS][j]/recip_fac,
 536                 dos[DOS_SOLID][j]/recip_fac,
 537                 dos[DOS_DIFF][j]/recip_fac);
 538     }
 539     xvgrclose(fp);
 540
 541     /* Finally analyze the results! */
 542     wCdiff = 0.5;
 543     wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
 544     wEdiff = 0.5;
 545     wAdiff = wEdiff-wSdiff;
 546     for (j = 0; (j < nframes/4); j++)
 547     {
 548         dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
 549                           dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
 550         dos[DOS_S][j]  = (dos[DOS_DIFF][j]*wSdiff +
 551                           dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
 552         dos[DOS_A][j]  = (dos[DOS_DIFF][j]*wAdiff +
 553                           dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
 554         dos[DOS_E][j]  = (dos[DOS_DIFF][j]*wEdiff +
 555                           dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
 556     }
 557     DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
 558     DiffCoeff = 1000*DiffCoeff/3.0;
 559     fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
 560             DiffCoeff);
 561     fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
 562             1000*DoS0/(12*tmass*beta));
 563
 564     cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
 565                                    nframes/4, &stddev);
 566     fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
 567
 568     /*
 569        S  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_S],NULL,
 570                                    nframes/4,&stddev);
 571        fprintf(fplog,"Entropy %g J/mol K\n",1000*S/Nmol);
 572        A  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_A],NULL,
 573                                    nframes/4,&stddev);
 574        fprintf(fplog,"Helmholtz energy %g kJ/mol\n",A/Nmol);
 575        E  = BOLTZ * evaluate_integral(nframes/4,nu,dos[DOS_E],NULL,
 576                                    nframes/4,&stddev);
 577        fprintf(fplog,"Internal energy %g kJ/mol\n",E/Nmol);
 578      */
 579     fprintf(fplog, "\nArrivederci!\n");
 580     gmx_fio_fclose(fplog);
 581
 582     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");
 583
 584     return 0;
 585 }