From: Mark Abraham <mark.j.abraham@gmail.com>
Date: Fri, 8 Jul 2016 13:41:00 +0000 (+0200)
Subject: Merge branch release-5-1 into release-2016
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=612274586adf48f4131541d11abf9e92837ac77b;p=alexxy%2Fgromacs.git

Merge branch release-5-1 into release-2016

Change-Id: I976ad584d1c8757c6464aae482cb5a4c245beacb
---

612274586adf48f4131541d11abf9e92837ac77b
diff --cc src/gromacs/gmxana/gmx_dos.cpp
index 9737716eb8,0000000000..88a3667a36
mode 100644,000000..100644
--- a/src/gromacs/gmxana/gmx_dos.cpp
+++ b/src/gromacs/gmxana/gmx_dos.cpp
@@@ -1,598 -1,0 +1,603 @@@
 +/*
 + * This file is part of the GROMACS molecular simulation package.
 + *
 + * Copyright (c) 2011,2012,2013,2014,2015,2016, by the GROMACS development team, led by
 + * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
 + * and including many others, as listed in the AUTHORS file in the
 + * top-level source directory and at http://www.gromacs.org.
 + *
 + * GROMACS is free software; you can redistribute it and/or
 + * modify it under the terms of the GNU Lesser General Public License
 + * as published by the Free Software Foundation; either version 2.1
 + * of the License, or (at your option) any later version.
 + *
 + * GROMACS is distributed in the hope that it will be useful,
 + * but WITHOUT ANY WARRANTY; without even the implied warranty of
 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 + * Lesser General Public License for more details.
 + *
 + * You should have received a copy of the GNU Lesser General Public
 + * License along with GROMACS; if not, see
 + * http://www.gnu.org/licenses, or write to the Free Software Foundation,
 + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
 + *
 + * If you want to redistribute modifications to GROMACS, please
 + * consider that scientific software is very special. Version
 + * control is crucial - bugs must be traceable. We will be happy to
 + * consider code for inclusion in the official distribution, but
 + * derived work must not be called official GROMACS. Details are found
 + * in the README & COPYING files - if they are missing, get the
 + * official version at http://www.gromacs.org.
 + *
 + * To help us fund GROMACS development, we humbly ask that you cite
 + * the research papers on the package. Check out http://www.gromacs.org.
 + */
 +#include "gmxpre.h"
 +
 +#include <cmath>
 +#include <cstdio>
 +#include <cstdlib>
 +#include <cstring>
 +
 +#include "gromacs/commandline/pargs.h"
 +#include "gromacs/commandline/viewit.h"
 +#include "gromacs/correlationfunctions/autocorr.h"
 +#include "gromacs/correlationfunctions/integrate.h"
 +#include "gromacs/fft/fft.h"
 +#include "gromacs/fileio/confio.h"
 +#include "gromacs/fileio/gmxfio.h"
 +#include "gromacs/fileio/trxio.h"
 +#include "gromacs/fileio/xvgr.h"
 +#include "gromacs/gmxana/gmx_ana.h"
 +#include "gromacs/math/functions.h"
 +#include "gromacs/math/units.h"
 +#include "gromacs/math/utilities.h"
 +#include "gromacs/math/vec.h"
 +#include "gromacs/topology/index.h"
 +#include "gromacs/topology/topology.h"
 +#include "gromacs/trajectory/trajectoryframe.h"
 +#include "gromacs/utility/arraysize.h"
 +#include "gromacs/utility/fatalerror.h"
 +#include "gromacs/utility/futil.h"
 +#include "gromacs/utility/pleasecite.h"
 +#include "gromacs/utility/smalloc.h"
 +
 +enum {
 +    VACF, MVACF, DOS, DOS_SOLID, DOS_DIFF, DOS_CP, DOS_S, DOS_A, DOS_E, DOS_NR
 +};
 +
 +static int calcMoleculesInIndexGroup(const t_block *mols, int natoms, const int *index, int nindex)
 +{
 +    int   i    = 0;
 +    int   mol  = 0;
 +    int   nMol = 0;
 +    int   j;
 +
 +    while (i < nindex)
 +    {
 +        while (index[i] > mols->index[mol])
 +        {
 +            mol++;
 +            if (mol >= mols->nr)
 +            {
 +                gmx_fatal(FARGS, "Atom index out of range: %d", index[i]+1);
 +            }
 +        }
 +        for (j = mols->index[mol]; j < mols->index[mol+1]; j++)
 +        {
 +            if (index[i] != j)
 +            {
 +                gmx_fatal(FARGS, "The index group does not consist of whole molecules");
 +            }
 +            i++;
 +            if (i > natoms)
 +            {
 +                gmx_fatal(FARGS, "Index contains atom numbers larger than the topology");
 +            }
 +        }
 +        nMol++;
 +    }
 +    return nMol;
 +}
 +
 +static double FD(double Delta, double f)
 +{
 +    return (2*std::pow(Delta, -4.5)*std::pow(f, 7.5) -
 +            6*std::pow(Delta, -3.0)*std::pow(f, 5.0) -
 +            std::pow(Delta, -1.5)*std::pow(f, 3.5) +
 +            6*std::pow(Delta, -1.5)*std::pow(f, 2.5) +
 +            2*f - 2);
 +}
 +
 +static double YYY(double f, double y)
 +{
 +    return (2*gmx::power3(y*f) - gmx::square(f)*y*(1+6*y) +
 +            (2+6*y)*f - 2);
 +}
 +
 +static double calc_compress(double y)
 +{
 +    if (y == 1)
 +    {
 +        return 0;
 +    }
 +    return ((1+y+gmx::square(y)-gmx::power3(y))/(gmx::power3(1-y)));
 +}
 +
 +static double bisector(double Delta, double tol,
 +                       double ff0, double ff1,
 +                       double ff(double, double))
 +{
 +    double fd, f, f0, f1;
 +    double tolmin = 1e-8;
 +
 +    f0 = ff0;
 +    f1 = ff1;
 +    if (tol < tolmin)
 +    {
 +        fprintf(stderr, "Unrealistic tolerance %g for bisector. Setting it to %g\n", tol, tolmin);
 +        tol = tolmin;
 +    }
 +
 +    do
 +    {
 +        f   = (f0+f1)*0.5;
 +        fd  = ff(Delta, f);
 +        if (fd < 0)
 +        {
 +            f0 = f;
 +        }
 +        else if (fd > 0)
 +        {
 +            f1 = f;
 +        }
 +        else
 +        {
 +            return f;
 +        }
 +    }
 +    while ((f1-f0) > tol);
 +
 +    return f;
 +}
 +
 +static double calc_fluidicity(double Delta, double tol)
 +{
 +    return bisector(Delta, tol, 0, 1, FD);
 +}
 +
 +static double calc_y(double f, double Delta, double toler)
 +{
 +    double y1, y2;
 +
 +    y1 = std::pow(f/Delta, 1.5);
 +    y2 = bisector(f, toler, 0, 10000, YYY);
 +    if (std::abs((y1-y2)/(y1+y2)) > 100*toler)
 +    {
 +        fprintf(stderr, "Inconsistency computing y: y1 = %f, y2 = %f, using y1.\n",
 +                y1, y2);
 +    }
 +
 +    return y1;
 +}
 +
 +static double calc_Shs(double f, double y)
 +{
 +    double fy  = f*y;
 +
 +    return BOLTZ*(std::log(calc_compress(fy)) + fy*(3*fy-4)/gmx::square(1-fy));
 +}
 +
 +static real wCsolid(real nu, real beta)
 +{
 +    real bhn = beta*PLANCK*nu;
 +    real ebn, koko;
 +
 +    if (bhn == 0)
 +    {
 +        return 1.0;
 +    }
 +    else
 +    {
 +        ebn  = std::exp(bhn);
 +        koko = gmx::square(1-ebn);
 +        return gmx::square(bhn)*ebn/koko;
 +    }
 +}
 +
 +static real wSsolid(real nu, real beta)
 +{
 +    real bhn = beta*PLANCK*nu;
 +
 +    if (bhn == 0)
 +    {
 +        return 1;
 +    }
 +    else
 +    {
 +        return bhn/std::expm1(bhn) - std::log1p(-std::exp(-bhn));
 +    }
 +}
 +
 +static real wAsolid(real nu, real beta)
 +{
 +    real bhn = beta*PLANCK*nu;
 +
 +    if (bhn == 0)
 +    {
 +        return 0;
 +    }
 +    else
 +    {
 +        return std::log((1-std::exp(-bhn))/(std::exp(-bhn/2))) - std::log(bhn);
 +    }
 +}
 +
 +static real wEsolid(real nu, real beta)
 +{
 +    real bhn = beta*PLANCK*nu;
 +
 +    if (bhn == 0)
 +    {
 +        return 1;
 +    }
 +    else
 +    {
 +        return bhn/2 + bhn/std::expm1(bhn)-1;
 +    }
 +}
 +
 +int gmx_dos(int argc, char *argv[])
 +{
 +    const char         *desc[] = {
 +        "[THISMODULE] computes the Density of States from a simulations.",
 +        "In order for this to be meaningful the velocities must be saved",
 +        "in the trajecotry with sufficiently high frequency such as to cover",
 +        "all vibrations. For flexible systems that would be around a few fs",
 +        "between saving. Properties based on the DoS are printed on the",
 +        "standard output."
 +        "Note that the density of states is calculated from the mass-weighted",
 +        "autocorrelation, and by default only from the square of the real",
 +        "component rather than absolute value. This means the shape can differ",
 +        "substantially from the plain vibrational power spectrum you can",
 +        "calculate with gmx velacc."
 +    };
 +    const char         *bugs[] = {
 +        "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)."
 +    };
 +    FILE               *fp, *fplog;
 +    t_topology          top;
 +    int                 ePBC = -1;
 +    t_trxframe          fr;
 +    matrix              box;
 +    int                 gnx;
 +    real                t0, t1;
 +    t_trxstatus        *status;
 +    int                 nV, nframes, n_alloc, i, j, fftcode, Nmol, Natom;
 +    double              rho, dt, Vsum, V, tmass, dostot, dos2;
 +    real              **c1, **dos, mi, beta, bfac, *nu, *tt, stddev, c1j;
 +    gmx_output_env_t   *oenv;
 +    gmx_fft_t           fft;
 +    double              cP, DiffCoeff, Delta, f, y, z, sigHS, Shs, Sig, DoS0, recip_fac;
 +    double              wCdiff, wSdiff, wAdiff, wEdiff;
 +    int                 grpNatoms;
 +    int                *index;
 +    char               *grpname;
 +    double              invNormalize;
 +    gmx_bool            normalizeAutocorrelation;
 +
-     static     gmx_bool bVerbose = TRUE, bAbsolute = FALSE, bNormalizeDos = FALSE;
-     static     gmx_bool bRecip   = FALSE;
-     static     real     Temp     = 298.15, toler = 1e-6;
++    static     gmx_bool bVerbose   = TRUE, bAbsolute = FALSE, bNormalizeDos = FALSE;
++    static     gmx_bool bRecip     = FALSE;
++    static     real     Temp       = 298.15, toler = 1e-6;
++    int                 min_frames = 100;
 +
 +    t_pargs             pa[]     = {
 +        { "-v", FALSE, etBOOL, {&bVerbose},
 +          "Be loud and noisy." },
 +        { "-recip", FALSE, etBOOL, {&bRecip},
 +          "Use cm^-1 on X-axis instead of 1/ps for DoS plots." },
 +        { "-abs", FALSE, etBOOL, {&bAbsolute},
 +          "Use the absolute value of the Fourier transform of the VACF as the Density of States. Default is to use the real component only" },
 +        { "-normdos", FALSE, etBOOL, {&bNormalizeDos},
 +          "Normalize the DoS such that it adds up to 3N. This should usually not be necessary." },
 +        { "-T", FALSE, etREAL, {&Temp},
 +          "Temperature in the simulation" },
 +        { "-toler", FALSE, etREAL, {&toler},
 +          "[HIDDEN]Tolerance when computing the fluidicity using bisection algorithm" }
 +    };
 +
 +    t_filenm            fnm[] = {
 +        { efTRN, "-f",    NULL,    ffREAD  },
 +        { efTPR, "-s",    NULL,    ffREAD  },
 +        { efNDX, NULL,    NULL,    ffOPTRD },
 +        { efXVG, "-vacf", "vacf",  ffWRITE },
 +        { efXVG, "-mvacf", "mvacf", ffWRITE },
 +        { efXVG, "-dos",  "dos",   ffWRITE },
 +        { efLOG, "-g",    "dos",   ffWRITE },
 +    };
 +#define NFILE asize(fnm)
 +    int                 npargs;
 +    t_pargs            *ppa;
 +    const char         *DoSlegend[] = {
 +        "DoS(v)", "DoS(v)[Solid]", "DoS(v)[Diff]"
 +    };
 +
 +    npargs = asize(pa);
 +    ppa    = add_acf_pargs(&npargs, pa);
 +    if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
 +                           NFILE, fnm, npargs, ppa, asize(desc), desc,
 +                           asize(bugs), bugs, &oenv))
 +    {
 +        return 0;
 +    }
 +
 +    beta = 1/(Temp*BOLTZ);
 +
 +    fplog = gmx_fio_fopen(ftp2fn(efLOG, NFILE, fnm), "w");
 +    fprintf(fplog, "Doing density of states analysis based on trajectory.\n");
 +    please_cite(fplog, "Pascal2011a");
 +    please_cite(fplog, "Caleman2011b");
 +
 +    read_tps_conf(ftp2fn(efTPR, NFILE, fnm), &top, &ePBC, NULL, NULL, box, TRUE);
 +
 +    /* Handle index groups */
 +    get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &grpNatoms, &index, &grpname);
 +
 +    V     = det(box);
 +    tmass = 0;
 +    for (i = 0; i < grpNatoms; i++)
 +    {
 +        tmass += top.atoms.atom[index[i]].m;
 +    }
 +
 +    Natom = grpNatoms;
 +    Nmol  = calcMoleculesInIndexGroup(&top.mols, top.atoms.nr, index, grpNatoms);
 +    gnx   = Natom*DIM;
 +
 +    /* Correlation stuff */
 +    snew(c1, gnx);
 +    for (i = 0; (i < gnx); i++)
 +    {
 +        c1[i] = NULL;
 +    }
 +
 +    read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_V);
 +    t0 = fr.time;
 +
 +    n_alloc = 0;
 +    nframes = 0;
 +    Vsum    = 0;
 +    nV      = 0;
 +    do
 +    {
 +        if (fr.bBox)
 +        {
 +            V      = det(fr.box);
 +            Vsum  += V;
 +            nV++;
 +        }
 +        if (nframes >= n_alloc)
 +        {
 +            n_alloc += 100;
 +            for (i = 0; i < gnx; i++)
 +            {
 +                srenew(c1[i], n_alloc);
 +            }
 +        }
 +        for (i = 0; i < gnx; i += DIM)
 +        {
 +            c1[i+XX][nframes] = fr.v[index[i/DIM]][XX];
 +            c1[i+YY][nframes] = fr.v[index[i/DIM]][YY];
 +            c1[i+ZZ][nframes] = fr.v[index[i/DIM]][ZZ];
 +        }
 +
 +        t1 = fr.time;
 +
 +        nframes++;
 +    }
 +    while (read_next_frame(oenv, status, &fr));
 +
 +    close_trj(status);
 +
++    if (nframes < min_frames)
++    {
++        gmx_fatal(FARGS, "You need at least %d frames in the trajectory and you only have %d.", min_frames, nframes);
++    }
 +    dt = (t1-t0)/(nframes-1);
 +    if (nV > 0)
 +    {
 +        V = Vsum/nV;
 +    }
 +    if (bVerbose)
 +    {
 +        printf("Going to do %d fourier transforms of length %d. Hang on.\n",
 +               gnx, nframes);
 +    }
 +    /* Unfortunately the -normalize program option for the autocorrelation
 +     * function calculation is added as a hack with a static variable in the
 +     * autocorrelation.c source. That would work if we called the normal
 +     * do_autocorr(), but this routine overrides that by directly calling
 +     * the low-level functionality. That unfortunately leads to ignoring the
 +     * default value for the option (which is to normalize).
 +     * Since the absolute value seems to be important for the subsequent
 +     * analysis below, we detect the value directly from the option, calculate
 +     * the autocorrelation without normalization, and then apply the
 +     * normalization just to the autocorrelation output
 +     * (or not, if the user asked for a non-normalized autocorrelation).
 +     */
 +    normalizeAutocorrelation = opt2parg_bool("-normalize", npargs, ppa);
 +
 +    /* Note that we always disable normalization here, regardless of user settings */
 +    low_do_autocorr(NULL, oenv, NULL, nframes, gnx, nframes, c1, dt, eacNormal, 0, FALSE,
 +                    FALSE, FALSE, -1, -1, 0);
 +    snew(dos, DOS_NR);
 +    for (j = 0; (j < DOS_NR); j++)
 +    {
 +        snew(dos[j], nframes+4);
 +    }
 +
 +    if (bVerbose)
 +    {
 +        printf("Going to merge the ACFs into the mass-weighted and plain ACF\n");
 +    }
 +    for (i = 0; (i < gnx); i += DIM)
 +    {
 +        mi = top.atoms.atom[index[i/DIM]].m;
 +        for (j = 0; (j < nframes/2); j++)
 +        {
 +            c1j            = (c1[i+XX][j] + c1[i+YY][j] + c1[i+ZZ][j]);
 +            dos[VACF][j]  += c1j/Natom;
 +            dos[MVACF][j] += mi*c1j;
 +        }
 +    }
 +
 +    fp = xvgropen(opt2fn("-vacf", NFILE, fnm), "Velocity autocorrelation function",
 +                  "Time (ps)", "C(t)", oenv);
 +    snew(tt, nframes/2);
 +
 +    invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0;
 +
 +    for (j = 0; (j < nframes/2); j++)
 +    {
 +        tt[j] = j*dt;
 +        fprintf(fp, "%10g  %10g\n", tt[j], dos[VACF][j] * invNormalize);
 +    }
 +    xvgrclose(fp);
 +
 +    fp = xvgropen(opt2fn("-mvacf", NFILE, fnm), "Mass-weighted velocity autocorrelation function",
 +                  "Time (ps)", "C(t)", oenv);
 +
 +    invNormalize = normalizeAutocorrelation ? 1.0/dos[VACF][0] : 1.0;
 +
 +    for (j = 0; (j < nframes/2); j++)
 +    {
 +        fprintf(fp, "%10g  %10g\n", tt[j], dos[MVACF][j] * invNormalize);
 +    }
 +    xvgrclose(fp);
 +
 +    if ((fftcode = gmx_fft_init_1d_real(&fft, nframes/2,
 +                                        GMX_FFT_FLAG_NONE)) != 0)
 +    {
 +        gmx_fatal(FARGS, "gmx_fft_init_1d_real returned %d", fftcode);
 +    }
 +    if ((fftcode = gmx_fft_1d_real(fft, GMX_FFT_REAL_TO_COMPLEX,
 +                                   (void *)dos[MVACF], (void *)dos[DOS])) != 0)
 +    {
 +        gmx_fatal(FARGS, "gmx_fft_1d_real returned %d", fftcode);
 +    }
 +
 +    /* First compute the DoS */
 +    /* Magic factor of 8 included now. */
 +    bfac = 8*dt*beta/2;
 +    dos2 = 0;
 +    snew(nu, nframes/4);
 +    for (j = 0; (j < nframes/4); j++)
 +    {
 +        nu[j] = 2*j/(t1-t0);
 +        dos2 += gmx::square(dos[DOS][2*j]) + gmx::square(dos[DOS][2*j+1]);
 +        if (bAbsolute)
 +        {
 +            dos[DOS][j] = bfac*std::hypot(dos[DOS][2*j], dos[DOS][2*j+1]);
 +        }
 +        else
 +        {
 +            dos[DOS][j] = bfac*dos[DOS][2*j];
 +        }
 +    }
 +    /* Normalize it */
 +    dostot = evaluate_integral(nframes/4, nu, dos[DOS], NULL, nframes/4, &stddev);
 +    if (bNormalizeDos)
 +    {
 +        for (j = 0; (j < nframes/4); j++)
 +        {
 +            dos[DOS][j] *= 3*Natom/dostot;
 +        }
 +    }
 +
 +    /* Now analyze it */
 +    DoS0 = dos[DOS][0];
 +
 +    /* Note this eqn. is incorrect in Pascal2011a! */
 +    Delta = ((2*DoS0/(9*Natom))*std::sqrt(M_PI*BOLTZ*Temp*Natom/tmass)*
 +             std::pow((Natom/V), 1.0/3.0)*std::pow(6.0/M_PI, 2.0/3.0));
 +    f     = calc_fluidicity(Delta, toler);
 +    y     = calc_y(f, Delta, toler);
 +    z     = calc_compress(y);
 +    Sig   = BOLTZ*(5.0/2.0+std::log(2*M_PI*BOLTZ*Temp/(gmx::square(PLANCK))*V/(f*Natom)));
 +    Shs   = Sig+calc_Shs(f, y);
 +    rho   = (tmass*AMU)/(V*NANO*NANO*NANO);
 +    sigHS = std::cbrt(6*y*V/(M_PI*Natom));
 +
 +    fprintf(fplog, "System = \"%s\"\n", *top.name);
 +    fprintf(fplog, "Nmol = %d\n", Nmol);
 +    fprintf(fplog, "Natom = %d\n", Natom);
 +    fprintf(fplog, "dt = %g ps\n", dt);
 +    fprintf(fplog, "tmass = %g amu\n", tmass);
 +    fprintf(fplog, "V = %g nm^3\n", V);
 +    fprintf(fplog, "rho = %g g/l\n", rho);
 +    fprintf(fplog, "T = %g K\n", Temp);
 +    fprintf(fplog, "beta = %g mol/kJ\n", beta);
 +
 +    fprintf(fplog, "\nDoS parameters\n");
 +    fprintf(fplog, "Delta = %g\n", Delta);
 +    fprintf(fplog, "fluidicity = %g\n", f);
 +    fprintf(fplog, "hard sphere packing fraction = %g\n", y);
 +    fprintf(fplog, "hard sphere compressibility = %g\n", z);
 +    fprintf(fplog, "ideal gas entropy = %g\n", Sig);
 +    fprintf(fplog, "hard sphere entropy = %g\n", Shs);
 +    fprintf(fplog, "sigma_HS = %g nm\n", sigHS);
 +    fprintf(fplog, "DoS0 = %g\n", DoS0);
 +    fprintf(fplog, "Dos2 = %g\n", dos2);
 +    fprintf(fplog, "DoSTot = %g\n", dostot);
 +
 +    /* Now compute solid (2) and diffusive (3) components */
 +    fp = xvgropen(opt2fn("-dos", NFILE, fnm), "Density of states",
 +                  bRecip ? "E (cm\\S-1\\N)" : "\\f{12}n\\f{4} (1/ps)",
 +                  "\\f{4}S(\\f{12}n\\f{4})", oenv);
 +    xvgr_legend(fp, asize(DoSlegend), DoSlegend, oenv);
 +    recip_fac = bRecip ? (1e7/SPEED_OF_LIGHT) : 1.0;
 +    for (j = 0; (j < nframes/4); j++)
 +    {
 +        dos[DOS_DIFF][j]  = DoS0/(1+gmx::square(DoS0*M_PI*nu[j]/(6*f*Natom)));
 +        dos[DOS_SOLID][j] = dos[DOS][j]-dos[DOS_DIFF][j];
 +        fprintf(fp, "%10g  %10g  %10g  %10g\n",
 +                recip_fac*nu[j],
 +                dos[DOS][j]/recip_fac,
 +                dos[DOS_SOLID][j]/recip_fac,
 +                dos[DOS_DIFF][j]/recip_fac);
 +    }
 +    xvgrclose(fp);
 +
 +    /* Finally analyze the results! */
 +    wCdiff = 0.5;
 +    wSdiff = Shs/(3*BOLTZ); /* Is this correct? */
 +    wEdiff = 0.5;
 +    wAdiff = wEdiff-wSdiff;
 +    for (j = 0; (j < nframes/4); j++)
 +    {
 +        dos[DOS_CP][j] = (dos[DOS_DIFF][j]*wCdiff +
 +                          dos[DOS_SOLID][j]*wCsolid(nu[j], beta));
 +        dos[DOS_S][j]  = (dos[DOS_DIFF][j]*wSdiff +
 +                          dos[DOS_SOLID][j]*wSsolid(nu[j], beta));
 +        dos[DOS_A][j]  = (dos[DOS_DIFF][j]*wAdiff +
 +                          dos[DOS_SOLID][j]*wAsolid(nu[j], beta));
 +        dos[DOS_E][j]  = (dos[DOS_DIFF][j]*wEdiff +
 +                          dos[DOS_SOLID][j]*wEsolid(nu[j], beta));
 +    }
 +    DiffCoeff = evaluate_integral(nframes/2, tt, dos[VACF], NULL, nframes/2, &stddev);
 +    DiffCoeff = 1000*DiffCoeff/3.0;
 +    fprintf(fplog, "Diffusion coefficient from VACF %g 10^-5 cm^2/s\n",
 +            DiffCoeff);
 +    fprintf(fplog, "Diffusion coefficient from DoS %g 10^-5 cm^2/s\n",
 +            1000*DoS0/(12*tmass*beta));
 +
 +    cP = BOLTZ * evaluate_integral(nframes/4, nu, dos[DOS_CP], NULL,
 +                                   nframes/4, &stddev);
 +    fprintf(fplog, "Heat capacity %g J/mol K\n", 1000*cP/Nmol);
 +    fprintf(fplog, "\nArrivederci!\n");
 +    gmx_fio_fclose(fplog);
 +
 +    do_view(oenv, ftp2fn(efXVG, NFILE, fnm), "-nxy");
 +
 +    return 0;
 +}