src/gromacs/gmxana/gmx_tcaf.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cmath>
  41 #include <cstdio>
  42 #include <cstring>
  43
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/commandline/viewit.h"
  46 #include "gromacs/correlationfunctions/autocorr.h"
  47 #include "gromacs/correlationfunctions/expfit.h"
  48 #include "gromacs/fileio/confio.h"
  49 #include "gromacs/fileio/trxio.h"
  50 #include "gromacs/fileio/xvgr.h"
  51 #include "gromacs/gmxana/gmx_ana.h"
  52 #include "gromacs/gmxana/gstat.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/math/utilities.h"
  55 #include "gromacs/math/vec.h"
  56 #include "gromacs/pbcutil/pbc.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/gmxassert.h"
  64 #include "gromacs/utility/smalloc.h"
  65
  66 #define NK 24
  67 #define NPK 4
  68
  69 #define NKC 6
  70 #define NKC0 4
  71 static const int kset_c[NKC + 1] = { 0, 3, 9, 13, 16, 19, NK };
  72
  73 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
  74 static rvec v0[NK] = { { 1, 0, 0 },  { 0, 1, 0 },  { 0, 0, 1 },  { 1, 1, 0 },  { 1, -1, 0 },
  75                        { 1, 0, 1 },  { 1, 0, -1 }, { 0, 1, 1 },  { 0, 1, -1 }, { 1, 1, 1 },
  76                        { 1, 1, -1 }, { 1, -1, 1 }, { -1, 1, 1 }, { 2, 0, 0 },  { 0, 2, 0 },
  77                        { 0, 0, 2 },  { 3, 0, 0 },  { 0, 3, 0 },  { 0, 0, 3 },  { 4, 0, 0 },
  78                        { 0, 4, 0 },  { 0, 0, 4 } };
  79 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
  80 static rvec v1[NK] = { { 0, 1, 0 },  { 0, 0, 1 },  { 1, 0, 0 },  { 0, 0, 1 }, { 0, 0, 1 },
  81                        { 0, 1, 0 },  { 0, 1, 0 },  { 1, 0, 0 },  { 1, 0, 0 }, { 1, -1, 0 },
  82                        { 1, -1, 0 }, { 1, 0, -1 }, { 0, 1, -1 }, { 0, 1, 0 }, { 0, 0, 1 },
  83                        { 1, 0, 0 },  { 0, 1, 0 },  { 0, 0, 1 },  { 1, 0, 0 }, { 0, 1, 0 },
  84                        { 0, 0, 1 },  { 1, 0, 0 } };
  85 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
  86 static rvec v2[NK] = { { 0, 0, 1 },  { 1, 0, 0 }, { 0, 1, 0 },  { 1, -1, 0 }, { 1, 1, 0 },
  87                        { 1, 0, -1 }, { 1, 0, 1 }, { 0, 1, -1 }, { 0, 1, 1 },  { 1, 1, -2 },
  88                        { 1, 1, 2 },  { 1, 2, 1 }, { 2, 1, 1 },  { 0, 0, 1 },  { 1, 0, 0 },
  89                        { 0, 1, 0 },  { 0, 0, 1 }, { 1, 0, 0 },  { 0, 1, 0 },  { 0, 0, 1 },
  90                        { 1, 0, 0 },  { 0, 1, 0 } };
  91
  92 static void process_tcaf(int                     nframes,
  93                          real                    dt,
  94                          int                     nkc,
  95                          real**                  tc,
  96                          rvec*                   kfac,
  97                          real                    rho,
  98                          real                    wt,
  99                          const char*             fn_trans,
 100                          const char*             fn_tca,
 101                          const char*             fn_tc,
 102                          const char*             fn_tcf,
 103                          const char*             fn_cub,
 104                          const char*             fn_vk,
 105                          const gmx_output_env_t* oenv)
 106 {
 107     FILE * fp, *fp_vk, *fp_cub = nullptr;
 108     int    nk, ntc;
 109     real **tcaf, **tcafc = nullptr, eta, *sig;
 110     int    i, j, k, kc;
 111     int    ncorr;
 112     double fitparms[3];
 113
 114     nk  = kset_c[nkc];
 115     ntc = nk * NPK;
 116
 117     if (fn_trans)
 118     {
 119         fp = xvgropen(fn_trans, "Transverse Current", "Time (ps)", "TC (nm/ps)", oenv);
 120         for (i = 0; i < nframes; i++)
 121         {
 122             fprintf(fp, "%g", i * dt);
 123             for (j = 0; j < ntc; j++)
 124             {
 125                 fprintf(fp, " %g", tc[j][i]);
 126             }
 127             fprintf(fp, "\n");
 128         }
 129         xvgrclose(fp);
 130         do_view(oenv, fn_trans, "-nxy");
 131     }
 132
 133     ncorr = (nframes + 1) / 2;
 134     if (ncorr > gmx::roundToInt(5 * wt / dt))
 135     {
 136         ncorr = gmx::roundToInt(5 * wt / dt) + 1;
 137     }
 138     snew(tcaf, nk);
 139     for (k = 0; k < nk; k++)
 140     {
 141         snew(tcaf[k], ncorr);
 142     }
 143     if (fn_cub)
 144     {
 145         snew(tcafc, nkc);
 146         for (k = 0; k < nkc; k++)
 147         {
 148             snew(tcafc[k], ncorr);
 149         }
 150     }
 151     snew(sig, ncorr);
 152     for (i = 0; i < ncorr; i++)
 153     {
 154         sig[i] = std::exp(0.5 * i * dt / wt);
 155     }
 156
 157     low_do_autocorr(fn_tca,
 158                     oenv,
 159                     "Transverse Current Autocorrelation Functions",
 160                     nframes,
 161                     ntc,
 162                     ncorr,
 163                     tc,
 164                     dt,
 165                     eacNormal,
 166                     1,
 167                     FALSE,
 168                     FALSE,
 169                     FALSE,
 170                     0,
 171                     0,
 172                     0);
 173     do_view(oenv, fn_tca, "-nxy");
 174
 175     fp = xvgropen(fn_tc, "Transverse Current Autocorrelation Functions", "Time (ps)", "TCAF", oenv);
 176     for (i = 0; i < ncorr; i++)
 177     {
 178         kc = 0;
 179         fprintf(fp, "%g", i * dt);
 180         for (k = 0; k < nk; k++)
 181         {
 182             for (j = 0; j < NPK; j++)
 183             {
 184                 tcaf[k][i] += tc[NPK * k + j][i];
 185             }
 186             if (fn_cub)
 187             {
 188                 for (j = 0; j < NPK; j++)
 189                 {
 190                     tcafc[kc][i] += tc[NPK * k + j][i];
 191                 }
 192             }
 193             if (i == 0)
 194             {
 195                 fprintf(fp, " %g", 1.0);
 196             }
 197             else
 198             {
 199                 tcaf[k][i] /= tcaf[k][0];
 200                 fprintf(fp, " %g", tcaf[k][i]);
 201             }
 202             if (k + 1 == kset_c[kc + 1])
 203             {
 204                 kc++;
 205             }
 206         }
 207         fprintf(fp, "\n");
 208     }
 209     xvgrclose(fp);
 210     do_view(oenv, fn_tc, "-nxy");
 211
 212     if (fn_cub)
 213     {
 214         fp_cub = xvgropen(fn_cub, "TCAFs and fits", "Time (ps)", "TCAF", oenv);
 215         for (kc = 0; kc < nkc; kc++)
 216         {
 217             fprintf(fp_cub, "%g %g\n", 0.0, 1.0);
 218             for (i = 1; i < ncorr; i++)
 219             {
 220                 tcafc[kc][i] /= tcafc[kc][0];
 221                 fprintf(fp_cub, "%g %g\n", i * dt, tcafc[kc][i]);
 222             }
 223             fprintf(fp_cub, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 224             tcafc[kc][0] = 1.0;
 225         }
 226     }
 227
 228     fp_vk = xvgropen(fn_vk, "Fits", "k (nm\\S-1\\N)", "\\8h\\4 (10\\S-3\\N kg m\\S-1\\N s\\S-1\\N)", oenv);
 229     if (output_env_get_print_xvgr_codes(oenv))
 230     {
 231         fprintf(fp_vk, "@    s0 symbol 2\n");
 232         fprintf(fp_vk, "@    s0 symbol color 1\n");
 233         fprintf(fp_vk, "@    s0 linestyle 0\n");
 234         if (fn_cub)
 235         {
 236             fprintf(fp_vk, "@    s1 symbol 3\n");
 237             fprintf(fp_vk, "@    s1 symbol color 2\n");
 238         }
 239     }
 240     fp = xvgropen(fn_tcf, "TCAF Fits", "Time (ps)", "", oenv);
 241     for (k = 0; k < nk; k++)
 242     {
 243         tcaf[k][0]  = 1.0;
 244         fitparms[0] = 1;
 245         fitparms[1] = 1;
 246         do_lmfit(ncorr, tcaf[k], sig, dt, nullptr, 0, ncorr * dt, oenv, bDebugMode(), effnVAC, fitparms, 0, nullptr);
 247         eta = 1000 * fitparms[1] * rho
 248               / (4 * fitparms[0] * gmx::c_pico * norm2(kfac[k]) / (gmx::c_nano * gmx::c_nano));
 249         fprintf(stdout, "k %6.3f  tau %6.3f  eta %8.5f 10^-3 kg/(m s)\n", norm(kfac[k]), fitparms[0], eta);
 250         fprintf(fp_vk, "%6.3f %g\n", norm(kfac[k]), eta);
 251         for (i = 0; i < ncorr; i++)
 252         {
 253             fprintf(fp, "%g %g\n", i * dt, fit_function(effnVAC, fitparms, i * dt));
 254         }
 255         fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 256     }
 257     xvgrclose(fp);
 258     do_view(oenv, fn_tcf, "-nxy");
 259
 260     if (fn_cub)
 261     {
 262         fprintf(stdout, "Averaged over k-vectors:\n");
 263         fprintf(fp_vk, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 264         for (k = 0; k < nkc; k++)
 265         {
 266             tcafc[k][0] = 1.0;
 267             fitparms[0] = 1;
 268             fitparms[1] = 1;
 269             do_lmfit(ncorr, tcafc[k], sig, dt, nullptr, 0, ncorr * dt, oenv, bDebugMode(), effnVAC, fitparms, 0, nullptr);
 270             eta = 1000 * fitparms[1] * rho
 271                   / (4 * fitparms[0] * gmx::c_pico * norm2(kfac[kset_c[k]]) / (gmx::c_nano * gmx::c_nano));
 272             fprintf(stdout,
 273                     "k %6.3f  tau %6.3f  Omega %6.3f  eta %8.5f 10^-3 kg/(m s)\n",
 274                     norm(kfac[kset_c[k]]),
 275                     fitparms[0],
 276                     fitparms[1],
 277                     eta);
 278             fprintf(fp_vk, "%6.3f %g\n", norm(kfac[kset_c[k]]), eta);
 279             for (i = 0; i < ncorr; i++)
 280             {
 281                 fprintf(fp_cub, "%g %g\n", i * dt, fit_function(effnVAC, fitparms, i * dt));
 282             }
 283             fprintf(fp_cub, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 284         }
 285         fprintf(fp_vk, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 286         xvgrclose(fp_cub);
 287         do_view(oenv, fn_cub, "-nxy");
 288     }
 289     xvgrclose(fp_vk);
 290     do_view(oenv, fn_vk, "-nxy");
 291 }
 292
 293
 294 int gmx_tcaf(int argc, char* argv[])
 295 {
 296     const char* desc[] = {
 297         "[THISMODULE] computes tranverse current autocorrelations.",
 298         "These are used to estimate the shear viscosity, [GRK]eta[grk].",
 299         "For details see: Palmer, Phys. Rev. E 49 (1994) pp 359-366.[PAR]",
 300         "Transverse currents are calculated using the",
 301         "k-vectors (1,0,0) and (2,0,0) each also in the [IT]y[it]- and [IT]z[it]-direction,",
 302         "(1,1,0) and (1,-1,0) each also in the 2 other planes (these vectors",
 303         "are not independent) and (1,1,1) and the 3 other box diagonals (also",
 304         "not independent). For each k-vector the sine and cosine are used, in",
 305         "combination with the velocity in 2 perpendicular directions. This gives",
 306         "a total of 16*2*2=64 transverse currents. One autocorrelation is",
 307         "calculated fitted for each k-vector, which gives 16 TCAFs. Each of",
 308         "these TCAFs is fitted to [MATH]f(t) = [EXP]-v[exp]([COSH]Wv[cosh] + 1/W ",
 309         "[SINH]Wv[sinh])[math],",
 310         "[MATH]v = -t/(2 [GRK]tau[grk])[math], [MATH]W = [SQRT]1 - 4 [GRK]tau[grk] ",
 311         "[GRK]eta[grk]/[GRK]rho[grk] k^2[sqrt][math], which gives 16 values of [GRK]tau[grk]",
 312         "and [GRK]eta[grk]. The fit weights decay exponentially with time constant [MATH]w[math] ",
 313         "(given with [TT]-wt[tt]) as [MATH][EXP]-t/w[exp][math], and the TCAF and",
 314         "fit are calculated up to time [MATH]5*w[math].",
 315         "The [GRK]eta[grk] values should be fitted to [MATH]1 - a [GRK]eta[grk](k) k^2[math], ",
 316         "from which one can estimate the shear viscosity at k=0.[PAR]",
 317         "When the box is cubic, one can use the option [TT]-oc[tt], which",
 318         "averages the TCAFs over all k-vectors with the same length.",
 319         "This results in more accurate TCAFs.",
 320         "Both the cubic TCAFs and fits are written to [TT]-oc[tt]",
 321         "The cubic [GRK]eta[grk] estimates are also written to [TT]-ov[tt].[PAR]",
 322         "With option [TT]-mol[tt], the transverse current is determined of",
 323         "molecules instead of atoms. In this case, the index group should",
 324         "consist of molecule numbers instead of atom numbers.[PAR]",
 325         "The k-dependent viscosities in the [TT]-ov[tt] file should be",
 326         "fitted to [MATH][GRK]eta[grk](k) = [GRK]eta[grk][SUB]0[sub] (1 - a k^2)[math] to obtain ",
 327         "the viscosity at",
 328         "infinite wavelength.[PAR]",
 329         "[BB]Note:[bb] make sure you write coordinates and velocities often enough.",
 330         "The initial, non-exponential, part of the autocorrelation function",
 331         "is very important for obtaining a good fit."
 332     };
 333
 334     static gmx_bool bMol = FALSE, bK34 = FALSE;
 335     static real     wt   = 5;
 336     t_pargs         pa[] = {
 337         { "-mol", FALSE, etBOOL, { &bMol }, "Calculate TCAF of molecules" },
 338         { "-k34", FALSE, etBOOL, { &bK34 }, "Also use k=(3,0,0) and k=(4,0,0)" },
 339         { "-wt", FALSE, etREAL, { &wt }, "Exponential decay time for the TCAF fit weights" }
 340     };
 341
 342     t_topology        top;
 343     PbcType           pbcType;
 344     t_trxframe        fr;
 345     matrix            box;
 346     gmx_bool          bTop;
 347     int               gnx;
 348     int *             index, *atndx = nullptr, at;
 349     char*             grpname;
 350     char              title[256];
 351     real              t0, t1, dt, m, mtot, sysmass, rho, sx, cx;
 352     t_trxstatus*      status;
 353     int               nframes, n_alloc, i, j, k, d;
 354     rvec              mv_mol, cm_mol, kfac[NK];
 355     int               nkc, nk, ntc;
 356     real**            tc;
 357     gmx_output_env_t* oenv;
 358
 359 #define NHISTO 360
 360
 361     t_filenm fnm[] = { { efTRN, "-f", nullptr, ffREAD },      { efTPS, nullptr, nullptr, ffOPTRD },
 362                        { efNDX, nullptr, nullptr, ffOPTRD },  { efXVG, "-ot", "transcur", ffOPTWR },
 363                        { efXVG, "-oa", "tcaf_all", ffWRITE }, { efXVG, "-o", "tcaf", ffWRITE },
 364                        { efXVG, "-of", "tcaf_fit", ffWRITE }, { efXVG, "-oc", "tcaf_cub", ffOPTWR },
 365                        { efXVG, "-ov", "visc_k", ffWRITE } };
 366 #define NFILE asize(fnm)
 367     int      npargs;
 368     t_pargs* ppa;
 369
 370     npargs = asize(pa);
 371     ppa    = add_acf_pargs(&npargs, pa);
 372
 373     if (!parse_common_args(
 374                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
 375     {
 376         sfree(ppa);
 377         return 0;
 378     }
 379
 380     bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, nullptr, nullptr, box, TRUE);
 381     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
 382
 383     if (bMol)
 384     {
 385         if (!bTop)
 386         {
 387             gmx_fatal(FARGS, "Need a topology to determine the molecules");
 388         }
 389         atndx = top.mols.index;
 390     }
 391
 392     if (bK34)
 393     {
 394         nkc = NKC;
 395     }
 396     else
 397     {
 398         nkc = NKC0;
 399     }
 400     nk = kset_c[nkc];
 401     GMX_ASSERT(nk >= 16, "Has to be over 16 because nkc is either NKC or NKC0.");
 402     ntc = nk * NPK;
 403
 404     sprintf(title, "Velocity Autocorrelation Function for %s", grpname);
 405
 406     sysmass = 0;
 407     for (i = 0; i < nk; i++)
 408     {
 409         if (iprod(v0[i], v1[i]) != 0)
 410         {
 411             gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
 412         }
 413         if (iprod(v0[i], v2[i]) != 0)
 414         {
 415             gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
 416         }
 417         if (iprod(v1[i], v2[i]) != 0)
 418         {
 419             gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
 420         }
 421         unitv(v1[i], v1[i]);
 422         unitv(v2[i], v2[i]);
 423     }
 424     snew(tc, ntc);
 425     for (i = 0; i < top.atoms.nr; i++)
 426     {
 427         sysmass += top.atoms.atom[i].m;
 428     }
 429
 430     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_X | TRX_NEED_V);
 431     t0 = fr.time;
 432
 433     n_alloc = 0;
 434     nframes = 0;
 435     rho     = 0;
 436
 437     do
 438     {
 439
 440         if (nframes >= n_alloc)
 441         {
 442             n_alloc += 100;
 443             for (i = 0; i < ntc; i++)
 444             {
 445                 srenew(tc[i], n_alloc);
 446             }
 447         }
 448
 449         rho += 1 / det(fr.box);
 450         for (k = 0; k < nk; k++)
 451         {
 452             for (d = 0; d < DIM; d++)
 453             {
 454                 kfac[k][d] = 2 * M_PI * v0[k][d] / fr.box[d][d];
 455             }
 456         }
 457         for (i = 0; i < ntc; i++)
 458         {
 459             tc[i][nframes] = 0;
 460         }
 461
 462         for (i = 0; i < gnx; i++)
 463         {
 464             if (bMol)
 465             {
 466                 clear_rvec(mv_mol);
 467                 clear_rvec(cm_mol);
 468                 mtot = 0;
 469                 for (j = 0; j < atndx[index[i] + 1] - atndx[index[i]]; j++)
 470                 {
 471                     at = atndx[index[i]] + j;
 472                     m  = top.atoms.atom[at].m;
 473                     mv_mol[XX] += m * fr.v[at][XX];
 474                     mv_mol[YY] += m * fr.v[at][YY];
 475                     mv_mol[ZZ] += m * fr.v[at][ZZ];
 476                     cm_mol[XX] += m * fr.x[at][XX];
 477                     cm_mol[YY] += m * fr.x[at][YY];
 478                     cm_mol[ZZ] += m * fr.x[at][ZZ];
 479                     mtot += m;
 480                 }
 481                 svmul(1.0 / mtot, cm_mol, cm_mol);
 482             }
 483             else
 484             {
 485                 svmul(top.atoms.atom[index[i]].m, fr.v[index[i]], mv_mol);
 486             }
 487
 488             if (!bMol)
 489             {
 490                 copy_rvec(fr.x[index[i]], cm_mol);
 491             }
 492             j = 0;
 493             for (k = 0; k < nk; k++)
 494             {
 495                 sx = std::sin(iprod(kfac[k], cm_mol));
 496                 cx = std::cos(iprod(kfac[k], cm_mol));
 497                 tc[j][nframes] += sx * iprod(v1[k], mv_mol);
 498                 j++;
 499                 tc[j][nframes] += cx * iprod(v1[k], mv_mol);
 500                 j++;
 501                 tc[j][nframes] += sx * iprod(v2[k], mv_mol);
 502                 j++;
 503                 tc[j][nframes] += cx * iprod(v2[k], mv_mol);
 504                 j++;
 505             }
 506         }
 507
 508         t1 = fr.time;
 509         nframes++;
 510     } while (read_next_frame(oenv, status, &fr));
 511     close_trx(status);
 512
 513     dt = (t1 - t0) / (nframes - 1);
 514
 515     rho *= sysmass / nframes * gmx::c_amu / (gmx::c_nano * gmx::c_nano * gmx::c_nano);
 516     fprintf(stdout, "Density = %g (kg/m^3)\n", rho);
 517     process_tcaf(nframes,
 518                  dt,
 519                  nkc,
 520                  tc,
 521                  kfac,
 522                  rho,
 523                  wt,
 524                  opt2fn_null("-ot", NFILE, fnm),
 525                  opt2fn("-oa", NFILE, fnm),
 526                  opt2fn("-o", NFILE, fnm),
 527                  opt2fn("-of", NFILE, fnm),
 528                  opt2fn_null("-oc", NFILE, fnm),
 529                  opt2fn("-ov", NFILE, fnm),
 530                  oenv);
 531
 532     return 0;
 533 }