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