src/gromacs/gmxana/gmx_dos.cpp

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