src/gromacs/gmxana/gmx_dos.cpp

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