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,2020,2021, 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/pbcutil/pbc.h"
  58 #include "gromacs/topology/index.h"
  59 #include "gromacs/topology/topology.h"
  60 #include "gromacs/trajectory/trajectoryframe.h"
  61 #include "gromacs/utility/arraysize.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/futil.h"
  64 #include "gromacs/utility/pleasecite.h"
  65 #include "gromacs/utility/smalloc.h"
  66
  67 enum
  68 {
  69     VACF,
  70     MVACF,
  71     DOS,
  72     DOS_SOLID,
  73     DOS_DIFF,
  74     DOS_CP,
  75     DOS_S,
  76     DOS_A,
  77     DOS_E,
  78     DOS_NR
  79 };
  80
  81 static int calcMoleculesInIndexGroup(const t_block* mols, int natoms, const int* index, int nindex)
  82 {
  83     int i    = 0;
  84     int mol  = 0;
  85     int nMol = 0;
  86     int j;
  87
  88     while (i < nindex)
  89     {
  90         while (index[i] > mols->index[mol])
  91         {
  92             mol++;
  93             if (mol >= mols->nr)
  94             {
  95                 gmx_fatal(FARGS, "Atom index out of range: %d", index[i] + 1);
  96             }
  97         }
  98         for (j = mols->index[mol]; j < mols->index[mol + 1]; j++)
  99         {
 100             if (index[i] != j)
 101             {
 102                 gmx_fatal(FARGS, "The index group does not consist of whole molecules");
 103             }
 104             i++;
 105             if (i > natoms)
 106             {
 107                 gmx_fatal(FARGS, "Index contains atom numbers larger than the topology");
 108             }
 109         }
 110         nMol++;
 111     }
 112     return nMol;
 113 }
 114
 115 static double FD(double Delta, double f)
 116 {
 117     return (2 * std::pow(Delta, -4.5) * std::pow(f, 7.5)
 118             - 6 * std::pow(Delta, -3.0) * std::pow(f, 5.0) - std::pow(Delta, -1.5) * std::pow(f, 3.5)
 119             + 6 * std::pow(Delta, -1.5) * std::pow(f, 2.5) + 2 * f - 2);
 120 }
 121
 122 static double YYY(double f, double y)
 123 {
 124     return (2 * gmx::power3(y * f) - gmx::square(f) * y * (1 + 6 * y) + (2 + 6 * y) * f - 2);
 125 }
 126
 127 static double calc_compress(double y)
 128 {
 129     if (y == 1)
 130     {
 131         return 0;
 132     }
 133     return ((1 + y + gmx::square(y) - gmx::power3(y)) / (gmx::power3(1 - y)));
 134 }
 135
 136 static double bisector(double Delta, double tol, double ff0, double ff1, double ff(double, double))
 137 {
 138     double fd, f, f0, f1;
 139     double tolmin = 1e-8;
 140
 141     f0 = ff0;
 142     f1 = ff1;
 143     if (tol < tolmin)
 144     {
 145         fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
 146         tol = tolmin;
 147     }
 148
 149     do
 150     {
 151         f  = (f0 + f1) * 0.5;
 152         fd = ff(Delta, f);
 153         if (fd < 0)
 154         {
 155             f0 = f;
 156         }
 157         else if (fd > 0)
 158         {
 159             f1 = f;
 160         }
 161         else
 162         {
 163             return f;
 164         }
 165     } while ((f1 - f0) > tol);
 166
 167     return f;
 168 }
 169
 170 static double calc_fluidicity(double Delta, double tol)
 171 {
 172     return bisector(Delta, tol, 0, 1, FD);
 173 }
 174
 175 static double calc_y(double f, double Delta, double toler)
 176 {
 177     double y1, y2;
 178
 179     y1 = std::pow(f / Delta, 1.5);
 180     y2 = bisector(f, toler, 0, 10000, YYY);
 181     if (std::abs((y1 - y2) / (y1 + y2)) > 100 * toler)
 182     {
 183         fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n", y1, y2);
 184     }
 185
 186     return y1;
 187 }
 188
 189 static double calc_Shs(double f, double y)
 190 {
 191     double fy = f * y;
 192
 193     return gmx::c_boltz * (std::log(calc_compress(fy)) + fy * (3 * fy - 4) / gmx::square(1 - fy));
 194 }
 195
 196 static real wCsolid(real nu, real beta)
 197 {
 198     real bhn = beta * gmx::c_planck * nu;
 199     real ebn, koko;
 200
 201     if (bhn == 0)
 202     {
 203         return 1.0;
 204     }
 205     else
 206     {
 207         ebn  = std::exp(bhn);
 208         koko = gmx::square(1 - ebn);
 209         return gmx::square(bhn) * ebn / koko;
 210     }
 211 }
 212
 213 static real wSsolid(real nu, real beta)
 214 {
 215     real bhn = beta * gmx::c_planck * nu;
 216
 217     if (bhn == 0)
 218     {
 219         return 1;
 220     }
 221     else
 222     {
 223         return bhn / std::expm1(bhn) - std::log1p(-std::exp(-bhn));
 224     }
 225 }
 226
 227 static real wAsolid(real nu, real beta)
 228 {
 229     real bhn = beta * gmx::c_planck * nu;
 230
 231     if (bhn == 0)
 232     {
 233         return 0;
 234     }
 235     else
 236     {
 237         return std::log((1 - std::exp(-bhn)) / (std::exp(-bhn / 2))) - std::log(bhn);
 238     }
 239 }
 240
 241 static real wEsolid(real nu, real beta)
 242 {
 243     real bhn = beta * gmx::c_planck * nu;
 244
 245     if (bhn == 0)
 246     {
 247         return 1;
 248     }
 249     else
 250     {
 251         return bhn / 2 + bhn / std::expm1(bhn) - 1;
 252     }
 253 }
 254
 255 int gmx_dos(int argc, char* argv[])
 256 {
 257     const char* desc[] = { "[THISMODULE] 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                            "Note that the density of states is calculated from the mass-weighted",
 264                            "autocorrelation, and by default only from the square of the real",
 265                            "component rather than absolute value. This means the shape can differ",
 266                            "substantially from the plain vibrational power spectrum you can",
 267                            "calculate with gmx velacc." };
 268     const char* bugs[] = {
 269         "This program needs a lot of memory: total usage equals the number of atoms times "
 270         "3 times number of frames times 4 (or 8 when run in double precision)."
 271     };
 272     FILE *            fp, *fplog;
 273     t_topology        top;
 274     PbcType           pbcType = PbcType::Unset;
 275     t_trxframe        fr;
 276     matrix            box;
 277     int               gnx;
 278     real              t0, t1;
 279     t_trxstatus*      status;
 280     int               nV, nframes, n_alloc, i, j, fftcode, Nmol, Natom;
 281     double            rho, dt, Vsum, V, tmass, dostot, dos2;
 282     real **           c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
 283     gmx_output_env_t* oenv;
 284     gmx_fft_t         fft;
 285     double            cP, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
 286     double            wCdiff, wSdiff, wAdiff, wEdiff;
 287     int               grpNatoms;
 288     int*              index;
 289     char*             grpname;
 290     double            invNormalize;
 291     gmx_bool          normalizeAutocorrelation;
 292
 293     static gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalizeDos = FALSE;
 294     static gmx_bool bRecip = FALSE;
 295     static real     Temp = 298.15, toler = 1e-6;
 296     int             min_frames = 100;
 297
 298     t_pargs pa[] = {
 299         { "-v", FALSE, etBOOL, { &bVerbose }, "Be loud and noisy." },
 300         { "-recip",
 301           FALSE,
 302           etBOOL,
 303           { &bRecip },
 304           "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
 305         { "-abs",
 306           FALSE,
 307           etBOOL,
 308           { &bAbsolute },
 309           "Use the absolute value of the Fourier transform of the VACF as the Density of States. "
 310           "Default is to use the real component only" },
 311         { "-normdos",
 312           FALSE,
 313           etBOOL,
 314           { &bNormalizeDos },
 315           "Normalize the DoS such that it adds up to 3N. This should usually not be necessary." },
 316         { "-T", FALSE, etREAL, { &Temp }, "Temperature in the simulation" },
 317         { "-toler",
 318           FALSE,
 319           etREAL,
 320           { &toler },
 321           "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" }
 322     };
 323
 324     t_filenm fnm[] = {
 325         { efTRN, "-f", nullptr, ffREAD },      { efTPR, "-s", nullptr, ffREAD },
 326         { efNDX, nullptr, nullptr, ffOPTRD },  { efXVG, "-vacf", "vacf", ffWRITE },
 327         { efXVG, "-mvacf", "mvacf", ffWRITE }, { efXVG, "-dos", "dos", ffWRITE },
 328         { efLOG, "-g", "dos", ffWRITE },
 329     };
 330 #define NFILE asize(fnm)
 331     int         npargs;
 332     t_pargs*    ppa;
 333     const char* DoSlegend[] = { "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]" };
 334
 335     npargs = asize(pa);
 336     ppa    = add_acf_pargs(&npargs, pa);
 337     if (!parse_common_args(
 338                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, asize(bugs), bugs, &oenv))
 339     {
 340         sfree(ppa);
 341         return 0;
 342     }
 343
 344     beta = 1 / (Temp * gmx::c_boltz);
 345
 346     fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
 347     fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
 348     please_cite(fplog, "Pascal2011a");
 349     please_cite(fplog, "Caleman2011b");
 350
 351     read_tps_conf(ftp2fn(efTPR, NFILE, fnm), &top, &pbcType, nullptr, nullptr, box, TRUE);
 352
 353     /* Handle index groups */
 354     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &grpNatoms, &index, &grpname);
 355
 356     V     = det(box);
 357     tmass = 0;
 358     for (i = 0; i < grpNatoms; i++)
 359     {
 360         tmass += top.atoms.atom[index[i]].m;
 361     }
 362
 363     Natom = grpNatoms;
 364     Nmol  = calcMoleculesInIndexGroup(&top.mols, top.atoms.nr, index, grpNatoms);
 365     gnx   = Natom * DIM;
 366
 367     /* Correlation stuff */
 368     snew(c1, gnx);
 369     for (i = 0; (i < gnx); i++)
 370     {
 371         c1[i] = nullptr;
 372     }
 373
 374     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
 375     t0 = fr.time;
 376
 377     n_alloc = 0;
 378     nframes = 0;
 379     Vsum    = 0;
 380     nV      = 0;
 381     do
 382     {
 383         if (fr.bBox)
 384         {
 385             V = det(fr.box);
 386             Vsum += V;
 387             nV++;
 388         }
 389         if (nframes >= n_alloc)
 390         {
 391             n_alloc += 100;
 392             for (i = 0; i < gnx; i++)
 393             {
 394                 srenew(c1[i], n_alloc);
 395             }
 396         }
 397         for (i = 0; i < gnx; i += DIM)
 398         {
 399             c1[i + XX][nframes] = fr.v[index[i / DIM]][XX];
 400             c1[i + YY][nframes] = fr.v[index[i / DIM]][YY];
 401             c1[i + ZZ][nframes] = fr.v[index[i / DIM]][ZZ];
 402         }
 403
 404         t1 = fr.time;
 405
 406         nframes++;
 407     } while (read_next_frame(oenv, status, &fr));
 408
 409     close_trx(status);
 410
 411     if (nframes < min_frames)
 412     {
 413         gmx_fatal(FARGS, "You need at least %d frames in the trajectory and you only have %d.", 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(
 440             nullptr, oenv, nullptr, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE, 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(
 463             opt2fn("-vacf", NFILE, fnm), "Velocity autocorrelation function", "Time (ps)", "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),
 476                   "Mass-weighted velocity autocorrelation function",
 477                   "Time (ps)",
 478                   "C(t)",
 479                   oenv);
 480
 481     invNormalize = normalizeAutocorrelation ? 1.0 / dos[VACF][0] : 1.0;
 482
 483     for (j = 0; (j < nframes / 2); j++)
 484     {
 485         fprintf(fp, "%10g  %10g\n", tt[j], dos[MVACF][j] * invNormalize);
 486     }
 487     xvgrclose(fp);
 488
 489     if ((fftcode = gmx_fft_init_1d_real(&fft, nframes / 2, GMX_FFT_FLAG_NONE)) != 0)
 490     {
 491         gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
 492     }
 493     if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX, dos[MVACF], dos[DOS])) != 0)
 494     {
 495         gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
 496     }
 497
 498     /* First compute the DoS */
 499     /* Magic factor of 8 included now. */
 500     bfac = 8 * dt * beta / 2;
 501     dos2 = 0;
 502     snew(nu, nframes / 4);
 503     for (j = 0; (j < nframes / 4); j++)
 504     {
 505         nu[j] = 2 * j / (t1 - t0);
 506         dos2 += gmx::square(dos[DOS][2 * j]) + gmx::square(dos[DOS][2 * j + 1]);
 507         if (bAbsolute)
 508         {
 509             dos[DOS][j] = bfac * std::hypot(dos[DOS][2 * j], dos[DOS][2 * j + 1]);
 510         }
 511         else
 512         {
 513             dos[DOS][j] = bfac * dos[DOS][2 * j];
 514         }
 515     }
 516     /* Normalize it */
 517     dostot = evaluate_integral(nframes / 4, nu, dos[DOS], nullptr, int{ nframes / 4 }, &stddev);
 518     if (bNormalizeDos)
 519     {
 520         for (j = 0; (j < nframes / 4); j++)
 521         {
 522             dos[DOS][j] *= 3 * Natom / dostot;
 523         }
 524     }
 525
 526     /* Now analyze it */
 527     DoS0 = dos[DOS][0];
 528
 529     /* Note this eqn. is incorrect in Pascal2011a! */
 530     Delta = ((2 * DoS0 / (9 * Natom)) * std::sqrt(M_PI * gmx::c_boltz * Temp * Natom / tmass)
 531              * std::pow((Natom / V), 1.0 / 3.0) * std::pow(6.0 / M_PI, 2.0 / 3.0));
 532     f     = calc_fluidicity(Delta, toler);
 533     y     = calc_y(f, Delta, toler);
 534     z     = calc_compress(y);
 535     Sig   = gmx::c_boltz
 536           * (5.0 / 2.0
 537              + std::log(2 * M_PI * gmx::c_boltz * Temp / (gmx::square(gmx::c_planck)) * V / (f * Natom)));
 538     Shs   = Sig + calc_Shs(f, y);
 539     rho   = (tmass * gmx::c_amu) / (V * gmx::c_nano * gmx::c_nano * gmx::c_nano);
 540     sigHS = std::cbrt(6 * y * V / (M_PI * Natom));
 541
 542     fprintf(fplog, "System = \"%s\"\n", *top.name);
 543     fprintf(fplog, "Nmol = %d\n", Nmol);
 544     fprintf(fplog, "Natom = %d\n", Natom);
 545     fprintf(fplog, "dt = %g ps\n", dt);
 546     fprintf(fplog, "tmass = %g amu\n", tmass);
 547     fprintf(fplog, "V = %g nm^3\n", V);
 548     fprintf(fplog, "rho = %g g/l\n", rho);
 549     fprintf(fplog, "T = %g K\n", Temp);
 550     fprintf(fplog, "beta = %g mol/kJ\n", beta);
 551
 552     fprintf(fplog, "\nDoS parameters\n");
 553     fprintf(fplog, "Delta = %g\n", Delta);
 554     fprintf(fplog, "fluidicity = %g\n", f);
 555     fprintf(fplog, "hard sphere packing fraction = %g\n", y);
 556     fprintf(fplog, "hard sphere compressibility = %g\n", z);
 557     fprintf(fplog, "ideal gas entropy = %g\n", Sig);
 558     fprintf(fplog, "hard sphere entropy = %g\n", Shs);
 559     fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
 560     fprintf(fplog, "DoS0 = %g\n", DoS0);
 561     fprintf(fplog, "Dos2 = %g\n", dos2);
 562     fprintf(fplog, "DoSTot = %g\n", dostot);
 563
 564     /* Now compute solid (2) and diffusive (3) components */
 565     fp = xvgropen(opt2fn("-dos", NFILE, fnm),
 566                   "Density of states",
 567                   bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
 568                   "\\f{4}S(\\f{12}n\\f{4})",
 569                   oenv);
 570     xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
 571     recip_fac = bRecip ? (1e7 / gmx::c_speedOfLight) : 1.0;
 572     for (j = 0; (j < nframes / 4); j++)
 573     {
 574         dos[DOS_DIFF][j]  = DoS0 / (1 + gmx::square(DoS0 * M_PI * nu[j] / (6 * f * Natom)));
 575         dos[DOS_SOLID][j] = dos[DOS][j] - dos[DOS_DIFF][j];
 576         fprintf(fp,
 577                 "%10g  %10g  %10g  %10g\n",
 578                 recip_fac * nu[j],
 579                 dos[DOS][j] / recip_fac,
 580                 dos[DOS_SOLID][j] / recip_fac,
 581                 dos[DOS_DIFF][j] / recip_fac);
 582     }
 583     xvgrclose(fp);
 584
 585     /* Finally analyze the results! */
 586     wCdiff = 0.5;
 587     wSdiff = Shs / (3 * gmx::c_boltz); /* Is this correct? */
 588     wEdiff = 0.5;
 589     wAdiff = wEdiff - wSdiff;
 590     for (j = 0; (j < nframes / 4); j++)
 591     {
 592         dos[DOS_CP][j] = (dos[DOS_DIFF][j] * wCdiff + dos[DOS_SOLID][j] * wCsolid(nu[j], beta));
 593         dos[DOS_S][j]  = (dos[DOS_DIFF][j] * wSdiff + dos[DOS_SOLID][j] * wSsolid(nu[j], beta));
 594         dos[DOS_A][j]  = (dos[DOS_DIFF][j] * wAdiff + dos[DOS_SOLID][j] * wAsolid(nu[j], beta));
 595         dos[DOS_E][j]  = (dos[DOS_DIFF][j] * wEdiff + dos[DOS_SOLID][j] * wEsolid(nu[j], beta));
 596     }
 597     DiffCoeff = evaluate_integral(nframes / 2, tt, dos[VACF], nullptr, nframes / 2., &stddev);
 598     DiffCoeff = 1000 * DiffCoeff / 3.0;
 599     fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n", DiffCoeff);
 600     fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n", 1000 * DoS0 / (12 * tmass * beta));
 601
 602     cP = gmx::c_boltz
 603          * evaluate_integral(nframes / 4, nu, dos[DOS_CP], nullptr, int{ nframes / 4 }, &stddev);
 604     fprintf(fplog, "Heat capacity %g J/mol K\n", 1000 * cP / Nmol);
 605     fprintf(fplog, "\nArrivederci!\n");
 606     gmx_fio_fclose(fplog);
 607
 608     do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");
 609
 610     return 0;
 611 }