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,2015, 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 #include "gmxpre.h"
  36
  37 #include <math.h>
  38 #include <stdio.h>
  39 #include <stdlib.h>
  40 #include <string.h>
  41
  42 #include "gromacs/commandline/pargs.h"
  43 #include "gromacs/correlationfunctions/autocorr.h"
  44 #include "gromacs/correlationfunctions/integrate.h"
  45 #include "gromacs/fft/fft.h"
  46 #include "gromacs/fileio/confio.h"
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/gmxana/gmx_ana.h"
  51 #include "gromacs/legacyheaders/copyrite.h"
  52 #include "gromacs/legacyheaders/macros.h"
  53 #include "gromacs/legacyheaders/txtdump.h"
  54 #include "gromacs/legacyheaders/typedefs.h"
  55 #include "gromacs/legacyheaders/viewit.h"
  56 #include "gromacs/math/units.h"
  57 #include "gromacs/math/utilities.h"
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/topology/index.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/futil.h"
  62 #include "gromacs/utility/smalloc.h"
  63 #include "gmx_ana.h"
  64
  65 enum {
  66     VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR
  67 };
  68
  69 static int calcMoleculesInIndexGroup(t_block *mols, int natoms, atom_id *index, int nindex)
  70 {
  71     int   i    = 0;
  72     int   mol  = 0;
  73     int   nMol = 0;
  74     int   j;
  75
  76     while (i < nindex)
  77     {
  78         while (index[i] > mols->index[mol])
  79         {
  80             mol++;
  81             if (mol >= mols->nr)
  82             {
  83                 gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
  84             }
  85         }
  86         for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
  87         {
  88             if (index[i] != j)
  89             {
  90                 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
  91             }
  92             i++;
  93             if (i == natoms)
  94             {
  95                 gmx_fatal(FARGS, "Index contains atom numbers larger than the topology");
  96             }
  97         }
  98         nMol++;
  99     }
 100     return nMol;
 101 }
 102
 103 static double FD(double Delta, double f)
 104 {
 105     return (2*pow(Delta, -4.5)*pow(f, 7.5) -
 106             6*pow(Delta, -3)*pow(f, 5) -
 107             pow(Delta, -1.5)*pow(f, 3.5) +
 108             6*pow(Delta, -1.5)*pow(f, 2.5) +
 109             2*f - 2);
 110 }
 111
 112 static double YYY(double f, double y)
 113 {
 114     return (2*pow(y*f, 3) - sqr(f)*y*(1+6*y) +
 115             (2+6*y)*f - 2);
 116 }
 117
 118 static double calc_compress(double y)
 119 {
 120     if (y == 1)
 121     {
 122         return 0;
 123     }
 124     return ((1+y+sqr(y)-pow(y, 3))/(pow(1-y, 3)));
 125 }
 126
 127 static double bisector(double Delta, double tol,
 128                        double ff0, double ff1,
 129                        double ff(double, double))
 130 {
 131     double fd0, fd, fd1, f, f0, f1;
 132     double tolmin = 1e-8;
 133
 134     f0 = ff0;
 135     f1 = ff1;
 136     if (tol < tolmin)
 137     {
 138         fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
 139         tol = tolmin;
 140     }
 141
 142     do
 143     {
 144         fd0 = ff(Delta, f0);
 145         fd1 = ff(Delta, f1);
 146         f   = (f0+f1)*0.5;
 147         fd  = ff(Delta, f);
 148         if (fd < 0)
 149         {
 150             f0 = f;
 151         }
 152         else if (fd > 0)
 153         {
 154             f1 = f;
 155         }
 156         else
 157         {
 158             return f;
 159         }
 160     }
 161     while ((f1-f0) > tol);
 162
 163     return f;
 164 }
 165
 166 static double calc_fluidicity(double Delta, double tol)
 167 {
 168     return bisector(Delta, tol, 0, 1, FD);
 169 }
 170
 171 static double calc_y(double f, double Delta, double toler)
 172 {
 173     double y1, y2;
 174
 175     y1 = pow(f/Delta, 1.5);
 176     y2 = bisector(f, toler, 0, 10000, YYY);
 177     if (fabs((y1-y2)/(y1+y2)) > 100*toler)
 178     {
 179         fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
 180                 y1, y2);
 181     }
 182
 183     return y1;
 184 }
 185
 186 static double calc_Shs(double f, double y)
 187 {
 188     double fy  = f*y;
 189
 190     return BOLTZ*(log(calc_compress(fy)) + fy*(3*fy-4)/sqr(1-fy));
 191 }
 192
 193 static real wCsolid(real nu, real beta)
 194 {
 195     real bhn = beta*PLANCK*nu;
 196     real ebn, koko;
 197
 198     if (bhn == 0)
 199     {
 200         return 1.0;
 201     }
 202     else
 203     {
 204         ebn  = exp(bhn);
 205         koko = sqr(1-ebn);
 206         return sqr(bhn)*ebn/koko;
 207     }
 208 }
 209
 210 static real wSsolid(real nu, real beta)
 211 {
 212     real bhn = beta*PLANCK*nu;
 213
 214     if (bhn == 0)
 215     {
 216         return 1;
 217     }
 218     else
 219     {
 220         return bhn/gmx_expm1(bhn) - gmx_log1p(-exp(-bhn));
 221     }
 222 }
 223
 224 static real wAsolid(real nu, real beta)
 225 {
 226     real bhn = beta*PLANCK*nu;
 227
 228     if (bhn == 0)
 229     {
 230         return 0;
 231     }
 232     else
 233     {
 234         return log((1-exp(-bhn))/(exp(-bhn/2))) - log(bhn);
 235     }
 236 }
 237
 238 static real wEsolid(real nu, real beta)
 239 {
 240     real bhn = beta*PLANCK*nu;
 241
 242     if (bhn == 0)
 243     {
 244         return 1;
 245     }
 246     else
 247     {
 248         return bhn/2 + bhn/gmx_expm1(bhn)-1;
 249     }
 250 }
 251
 252 int gmx_dos(int argc, char *argv[])
 253 {
 254     const char         *desc[] = {
 255         "[THISMODULE] computes the Density of States from a simulations.",
 256         "In order for this to be meaningful the velocities must be saved",
 257         "in the trajecotry with sufficiently high frequency such as to cover",
 258         "all vibrations. For flexible systems that would be around a few fs",
 259         "between saving. Properties based on the DoS are printed on the",
 260         "standard output."
 261         "Note that the density of states is calculated from the mass-weighted",
 262         "autocorrelation, and by default only from the square of the real",
 263         "component rather than absolute value. This means the shape can differ",
 264         "substantially from the plain vibrational power spectrum you can",
 265         "calculate with gmx velacc."
 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     int                 grpNatoms;
 287     atom_id            *index;
 288     char               *grpname;
 289     double              invNormalize;
 290     gmx_bool            normalizeAutocorrelation;
 291
 292     static     gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalizeDos = FALSE;
 293     static     gmx_bool bRecip   = FALSE;
 294     static     real     Temp     = 298.15, toler = 1e-6;
 295
 296     t_pargs             pa[]     = {
 297         { "-v", FALSE, etBOOL, {&bVerbose},
 298           "Be loud and noisy." },
 299         { "-recip", FALSE, etBOOL, {&bRecip},
 300           "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
 301         { "-abs", FALSE, etBOOL, {&bAbsolute},
 302           "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
 303         { "-normdos", FALSE, etBOOL, {&bNormalizeDos},
 304           "Normalize the DoS such that it adds up to 3N. This should usually not be necessary." },
 305         { "-T", FALSE, etREAL, {&Temp},
 306           "Temperature in the simulation" },
 307         { "-toler", FALSE, etREAL, {&toler},
 308           "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" }
 309     };
 310
 311     t_filenm            fnm[] = {
 312         { efTRN, "-f",    NULL,    ffREAD  },
 313         { efTPR, "-s",    NULL,    ffREAD  },
 314         { efNDX, NULL,    NULL,    ffOPTRD },
 315         { efXVG, "-vacf", "vacf",  ffWRITE },
 316         { efXVG, "-mvacf", "mvacf", ffWRITE },
 317         { efXVG, "-dos",  "dos",   ffWRITE },
 318         { efLOG, "-g",    "dos",   ffWRITE },
 319     };
 320 #define NFILE asize(fnm)
 321     int                 npargs;
 322     t_pargs            *ppa;
 323     const char         *DoSlegend[] = {
 324         "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
 325     };
 326
 327     npargs = asize(pa);
 328     ppa    = add_acf_pargs(&npargs, pa);
 329     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 330                            NFILE, fnm, npargs, ppa, asize(desc), desc,
 331                            asize(bugs), bugs, &oenv))
 332     {
 333         return 0;
 334     }
 335
 336     beta = 1/(Temp*BOLTZ);
 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(efTPR, NFILE, fnm), title, &top, &ePBC, NULL, NULL, box, TRUE);
 344
 345     /* Handle index groups */
 346     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &grpNatoms, &index, &grpname);
 347
 348     V     = det(box);
 349     tmass = 0;
 350     for (i = 0; i < grpNatoms; i++)
 351     {
 352         tmass += top.atoms.atom[index[i]].m;
 353     }
 354
 355     Natom = grpNatoms;
 356     Nmol  = calcMoleculesInIndexGroup(&top.mols, top.atoms.nr, index, grpNatoms);
 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[index[i/DIM]][XX];
 393             c1[i+YY][nframes] = fr.v[index[i/DIM]][YY];
 394             c1[i+ZZ][nframes] = fr.v[index[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     /* Unfortunately the -normalize program option for the autocorrelation
 416      * function calculation is added as a hack with a static variable in the
 417      * autocorrelation.c source. That would work if we called the normal
 418      * do_autocorr(), but this routine overrides that by directly calling
 419      * the low-level functionality. That unfortunately leads to ignoring the
 420      * default value for the option (which is to normalize).
 421      * Since the absolute value seems to be important for the subsequent
 422      * analysis below, we detect the value directly from the option, calculate
 423      * the autocorrelation without normalization, and then apply the
 424      * normalization just to the autocorrelation output
 425      * (or not, if the user asked for a non-normalized autocorrelation).
 426      */
 427     normalizeAutocorrelation = opt2parg_bool("-normalize", npargs, ppa);
 428
 429     /* Note that we always disable normalization here, regardless of user settings */
 430     low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
 431                     FALSE, FALSE, -1, -1, 0);
 432     snew(dos, DOS_NR);
 433     for (j = 0; (j < DOS_NR); j++)
 434     {
 435         snew(dos[j], nframes+4);
 436     }
 437
 438     if (bVerbose)
 439     {
 440         printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
 441     }
 442     for (i = 0; (i < gnx); i += DIM)
 443     {
 444         mi = top.atoms.atom[index[i/DIM]].m;
 445         for (j = 0; (j < nframes/2); j++)
 446         {
 447             c1j            = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
 448             dos[VACF][j]  += c1j/Natom;
 449             dos[MVACF][j] += mi*c1j;
 450         }
 451     }
 452
 453     fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity autocorrelation function",
 454                   "Time (ps)", "C(t)", oenv);
 455     snew(tt, nframes/2);
 456
 457     invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0;
 458
 459     for (j = 0; (j < nframes/2); j++)
 460     {
 461         tt[j] = j*dt;
 462         fprintf(fp, "%10g  %10g\n", tt[j], dos[VACF][j] * invNormalize);
 463     }
 464     xvgrclose(fp);
 465
 466     fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity autocorrelation function",
 467                   "Time (ps)", "C(t)", oenv);
 468
 469     invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0;
 470
 471     for (j = 0; (j < nframes/2); j++)
 472     {
 473         fprintf(fp, "%10g  %10g\n", tt[j], dos[MVACF][j] * invNormalize);
 474     }
 475     xvgrclose(fp);
 476
 477     if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
 478                                         GMX_FFT_FLAG_NONE)) != 0)
 479     {
 480         gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
 481     }
 482     if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
 483                                    (void *)dos[MVACF], (void *)dos[DOS])) != 0)
 484     {
 485         gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
 486     }
 487
 488     /* First compute the DoS */
 489     /* Magic factor of 8 included now. */
 490     bfac = 8*dt*beta/2;
 491     dos2 = 0;
 492     snew(nu, nframes/4);
 493     for (j = 0; (j < nframes/4); j++)
 494     {
 495         nu[j] = 2*j/(t1-t0);
 496         dos2 += sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]);
 497         if (bAbsolute)
 498         {
 499             dos[DOS][j] = bfac*sqrt(sqr(dos[DOS][2*j]) + sqr(dos[DOS][2*j+1]));
 500         }
 501         else
 502         {
 503             dos[DOS][j] = bfac*dos[DOS][2*j];
 504         }
 505     }
 506     /* Normalize it */
 507     dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
 508     if (bNormalizeDos)
 509     {
 510         for (j = 0; (j < nframes/4); j++)
 511         {
 512             dos[DOS][j] *= 3*Natom/dostot;
 513         }
 514     }
 515
 516     /* Now analyze it */
 517     DoS0 = dos[DOS][0];
 518
 519     /* Note this eqn. is incorrect in Pascal2011a! */
 520     Delta = ((2*DoS0/(9*Natom))*sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
 521              pow((Natom/V), 1.0/3.0)*pow(6/M_PI, 2.0/3.0));
 522     f     = calc_fluidicity(Delta, toler);
 523     y     = calc_y(f, Delta, toler);
 524     z     = calc_compress(y);
 525     Sig   = BOLTZ*(5.0/2.0+log(2*M_PI*BOLTZ*Temp/(sqr(PLANCK))*V/(f*Natom)));
 526     Shs   = Sig+calc_Shs(f, y);
 527     rho   = (tmass*AMU)/(V*NANO*NANO*NANO);
 528     sigHS = pow(6*y*V/(M_PI*Natom), 1.0/3.0);
 529
 530     fprintf(fplog, "System = \"%s\"\n", title);
 531     fprintf(fplog, "Nmol = %d\n", Nmol);
 532     fprintf(fplog, "Natom = %d\n", Natom);
 533     fprintf(fplog, "dt = %g ps\n", dt);
 534     fprintf(fplog, "tmass = %g amu\n", tmass);
 535     fprintf(fplog, "V = %g nm^3\n", V);
 536     fprintf(fplog, "rho = %g g/l\n", rho);
 537     fprintf(fplog, "T = %g K\n", Temp);
 538     fprintf(fplog, "beta = %g mol/kJ\n", beta);
 539
 540     fprintf(fplog, "\nDoS parameters\n");
 541     fprintf(fplog, "Delta = %g\n", Delta);
 542     fprintf(fplog, "fluidicity = %g\n", f);
 543     fprintf(fplog, "hard sphere packing fraction = %g\n", y);
 544     fprintf(fplog, "hard sphere compressibility = %g\n", z);
 545     fprintf(fplog, "ideal gas entropy = %g\n", Sig);
 546     fprintf(fplog, "hard sphere entropy = %g\n", Shs);
 547     fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
 548     fprintf(fplog, "DoS0 = %g\n", DoS0);
 549     fprintf(fplog, "Dos2 = %g\n", dos2);
 550     fprintf(fplog, "DoSTot = %g\n", dostot);
 551
 552     /* Now compute solid (2) and diffusive (3) components */
 553     fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
 554                   bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
 555                   "\\f{4}S(\\f{12}n\\f{4})", oenv);
 556     xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
 557     recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
 558     for (j = 0; (j < nframes/4); j++)
 559     {
 560         dos[DOS_DIFF][j]  = DoS0/(1+sqr(DoS0*M_PI*nu[j]/(6*f*Natom)));
 561         dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
 562         fprintf(fp, "%10g  %10g  %10g  %10g\n",
 563                 recip_fac*nu[j],
 564                 dos[DOS][j]/recip_fac,
 565                 dos[DOS_SOLID][j]/recip_fac,
 566                 dos[DOS_DIFF][j]/recip_fac);
 567     }
 568     xvgrclose(fp);
 569
 570     /* Finally analyze the results! */
 571     wCdiff = 0.5;
 572     wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
 573     wEdiff = 0.5;
 574     wAdiff = wEdiff-wSdiff;
 575     for (j = 0; (j < nframes/4); j++)
 576     {
 577         dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
 578                           dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
 579         dos[DOS_S][j]  = (dos[DOS_DIFF][j]*wSdiff +
 580                           dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
 581         dos[DOS_A][j]  = (dos[DOS_DIFF][j]*wAdiff +
 582                           dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
 583         dos[DOS_E][j]  = (dos[DOS_DIFF][j]*wEdiff +
 584                           dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
 585     }
 586     DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
 587     DiffCoeff = 1000*DiffCoeff/3.0;
 588     fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
 589             DiffCoeff);
 590     fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
 591             1000*DoS0/(12*tmass*beta));
 592
 593     cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
 594                                    nframes/4, &stddev);
 595     fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
 596     fprintf(fplog, "\nArrivederci!\n");
 597     gmx_fio_fclose(fplog);
 598
 599     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");
 600
 601     return 0;
 602 }