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