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