a0d0ce2ae6f989cba0cb8c49572e4c50f347a592
[alexxy/gromacs.git] / 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, 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 static rvec v0[NK] = { { 1, 0, 0 },  { 0, 1, 0 },  { 0, 0, 1 },  { 1, 1, 0 },  { 1, -1, 0 },
74                        { 1, 0, 1 },  { 1, 0, -1 }, { 0, 1, 1 },  { 0, 1, -1 }, { 1, 1, 1 },
75                        { 1, 1, -1 }, { 1, -1, 1 }, { -1, 1, 1 }, { 2, 0, 0 },  { 0, 2, 0 },
76                        { 0, 0, 2 },  { 3, 0, 0 },  { 0, 3, 0 },  { 0, 0, 3 },  { 4, 0, 0 },
77                        { 0, 4, 0 },  { 0, 0, 4 } };
78 static rvec v1[NK] = { { 0, 1, 0 },  { 0, 0, 1 },  { 1, 0, 0 },  { 0, 0, 1 }, { 0, 0, 1 },
79                        { 0, 1, 0 },  { 0, 1, 0 },  { 1, 0, 0 },  { 1, 0, 0 }, { 1, -1, 0 },
80                        { 1, -1, 0 }, { 1, 0, -1 }, { 0, 1, -1 }, { 0, 1, 0 }, { 0, 0, 1 },
81                        { 1, 0, 0 },  { 0, 1, 0 },  { 0, 0, 1 },  { 1, 0, 0 }, { 0, 1, 0 },
82                        { 0, 0, 1 },  { 1, 0, 0 } };
83 static rvec v2[NK] = { { 0, 0, 1 },  { 1, 0, 0 }, { 0, 1, 0 },  { 1, -1, 0 }, { 1, 1, 0 },
84                        { 1, 0, -1 }, { 1, 0, 1 }, { 0, 1, -1 }, { 0, 1, 1 },  { 1, 1, -2 },
85                        { 1, 1, 2 },  { 1, 2, 1 }, { 2, 1, 1 },  { 0, 0, 1 },  { 1, 0, 0 },
86                        { 0, 1, 0 },  { 0, 0, 1 }, { 1, 0, 0 },  { 0, 1, 0 },  { 0, 0, 1 },
87                        { 1, 0, 0 },  { 0, 1, 0 } };
88
89 static void process_tcaf(int                     nframes,
90                          real                    dt,
91                          int                     nkc,
92                          real**                  tc,
93                          rvec*                   kfac,
94                          real                    rho,
95                          real                    wt,
96                          const char*             fn_trans,
97                          const char*             fn_tca,
98                          const char*             fn_tc,
99                          const char*             fn_tcf,
100                          const char*             fn_cub,
101                          const char*             fn_vk,
102                          const gmx_output_env_t* oenv)
103 {
104     FILE * fp, *fp_vk, *fp_cub = nullptr;
105     int    nk, ntc;
106     real **tcaf, **tcafc = nullptr, eta, *sig;
107     int    i, j, k, kc;
108     int    ncorr;
109     double fitparms[3];
110
111     nk  = kset_c[nkc];
112     ntc = nk * NPK;
113
114     if (fn_trans)
115     {
116         fp = xvgropen(fn_trans, "Transverse Current", "Time (ps)", "TC (nm/ps)", oenv);
117         for (i = 0; i < nframes; i++)
118         {
119             fprintf(fp, "%g", i * dt);
120             for (j = 0; j < ntc; j++)
121             {
122                 fprintf(fp, " %g", tc[j][i]);
123             }
124             fprintf(fp, "\n");
125         }
126         xvgrclose(fp);
127         do_view(oenv, fn_trans, "-nxy");
128     }
129
130     ncorr = (nframes + 1) / 2;
131     if (ncorr > gmx::roundToInt(5 * wt / dt))
132     {
133         ncorr = gmx::roundToInt(5 * wt / dt) + 1;
134     }
135     snew(tcaf, nk);
136     for (k = 0; k < nk; k++)
137     {
138         snew(tcaf[k], ncorr);
139     }
140     if (fn_cub)
141     {
142         snew(tcafc, nkc);
143         for (k = 0; k < nkc; k++)
144         {
145             snew(tcafc[k], ncorr);
146         }
147     }
148     snew(sig, ncorr);
149     for (i = 0; i < ncorr; i++)
150     {
151         sig[i] = std::exp(0.5 * i * dt / wt);
152     }
153
154     low_do_autocorr(fn_tca,
155                     oenv,
156                     "Transverse Current Autocorrelation Functions",
157                     nframes,
158                     ntc,
159                     ncorr,
160                     tc,
161                     dt,
162                     eacNormal,
163                     1,
164                     FALSE,
165                     FALSE,
166                     FALSE,
167                     0,
168                     0,
169                     0);
170     do_view(oenv, fn_tca, "-nxy");
171
172     fp = xvgropen(fn_tc, "Transverse Current Autocorrelation Functions", "Time (ps)", "TCAF", oenv);
173     for (i = 0; i < ncorr; i++)
174     {
175         kc = 0;
176         fprintf(fp, "%g", i * dt);
177         for (k = 0; k < nk; k++)
178         {
179             for (j = 0; j < NPK; j++)
180             {
181                 tcaf[k][i] += tc[NPK * k + j][i];
182             }
183             if (fn_cub)
184             {
185                 for (j = 0; j < NPK; j++)
186                 {
187                     tcafc[kc][i] += tc[NPK * k + j][i];
188                 }
189             }
190             if (i == 0)
191             {
192                 fprintf(fp, " %g", 1.0);
193             }
194             else
195             {
196                 tcaf[k][i] /= tcaf[k][0];
197                 fprintf(fp, " %g", tcaf[k][i]);
198             }
199             if (k + 1 == kset_c[kc + 1])
200             {
201                 kc++;
202             }
203         }
204         fprintf(fp, "\n");
205     }
206     xvgrclose(fp);
207     do_view(oenv, fn_tc, "-nxy");
208
209     if (fn_cub)
210     {
211         fp_cub = xvgropen(fn_cub, "TCAFs and fits", "Time (ps)", "TCAF", oenv);
212         for (kc = 0; kc < nkc; kc++)
213         {
214             fprintf(fp_cub, "%g %g\n", 0.0, 1.0);
215             for (i = 1; i < ncorr; i++)
216             {
217                 tcafc[kc][i] /= tcafc[kc][0];
218                 fprintf(fp_cub, "%g %g\n", i * dt, tcafc[kc][i]);
219             }
220             fprintf(fp_cub, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
221             tcafc[kc][0] = 1.0;
222         }
223     }
224
225     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);
226     if (output_env_get_print_xvgr_codes(oenv))
227     {
228         fprintf(fp_vk, "@    s0 symbol 2\n");
229         fprintf(fp_vk, "@    s0 symbol color 1\n");
230         fprintf(fp_vk, "@    s0 linestyle 0\n");
231         if (fn_cub)
232         {
233             fprintf(fp_vk, "@    s1 symbol 3\n");
234             fprintf(fp_vk, "@    s1 symbol color 2\n");
235         }
236     }
237     fp = xvgropen(fn_tcf, "TCAF Fits", "Time (ps)", "", oenv);
238     for (k = 0; k < nk; k++)
239     {
240         tcaf[k][0]  = 1.0;
241         fitparms[0] = 1;
242         fitparms[1] = 1;
243         do_lmfit(ncorr, tcaf[k], sig, dt, nullptr, 0, ncorr * dt, oenv, bDebugMode(), effnVAC, fitparms, 0, nullptr);
244         eta = 1000 * fitparms[1] * rho / (4 * fitparms[0] * PICO * norm2(kfac[k]) / (NANO * NANO));
245         fprintf(stdout, "k %6.3f  tau %6.3f  eta %8.5f 10^-3 kg/(m s)\n", norm(kfac[k]), fitparms[0], eta);
246         fprintf(fp_vk, "%6.3f %g\n", norm(kfac[k]), eta);
247         for (i = 0; i < ncorr; i++)
248         {
249             fprintf(fp, "%g %g\n", i * dt, fit_function(effnVAC, fitparms, i * dt));
250         }
251         fprintf(fp, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
252     }
253     xvgrclose(fp);
254     do_view(oenv, fn_tcf, "-nxy");
255
256     if (fn_cub)
257     {
258         fprintf(stdout, "Averaged over k-vectors:\n");
259         fprintf(fp_vk, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
260         for (k = 0; k < nkc; k++)
261         {
262             tcafc[k][0] = 1.0;
263             fitparms[0] = 1;
264             fitparms[1] = 1;
265             do_lmfit(ncorr, tcafc[k], sig, dt, nullptr, 0, ncorr * dt, oenv, bDebugMode(), effnVAC, fitparms, 0, nullptr);
266             eta = 1000 * fitparms[1] * rho
267                   / (4 * fitparms[0] * PICO * norm2(kfac[kset_c[k]]) / (NANO * NANO));
268             fprintf(stdout,
269                     "k %6.3f  tau %6.3f  Omega %6.3f  eta %8.5f 10^-3 kg/(m s)\n",
270                     norm(kfac[kset_c[k]]),
271                     fitparms[0],
272                     fitparms[1],
273                     eta);
274             fprintf(fp_vk, "%6.3f %g\n", norm(kfac[kset_c[k]]), eta);
275             for (i = 0; i < ncorr; i++)
276             {
277                 fprintf(fp_cub, "%g %g\n", i * dt, fit_function(effnVAC, fitparms, i * dt));
278             }
279             fprintf(fp_cub, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
280         }
281         fprintf(fp_vk, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
282         xvgrclose(fp_cub);
283         do_view(oenv, fn_cub, "-nxy");
284     }
285     xvgrclose(fp_vk);
286     do_view(oenv, fn_vk, "-nxy");
287 }
288
289
290 int gmx_tcaf(int argc, char* argv[])
291 {
292     const char* desc[] = {
293         "[THISMODULE] computes tranverse current autocorrelations.",
294         "These are used to estimate the shear viscosity, [GRK]eta[grk].",
295         "For details see: Palmer, Phys. Rev. E 49 (1994) pp 359-366.[PAR]",
296         "Transverse currents are calculated using the",
297         "k-vectors (1,0,0) and (2,0,0) each also in the [IT]y[it]- and [IT]z[it]-direction,",
298         "(1,1,0) and (1,-1,0) each also in the 2 other planes (these vectors",
299         "are not independent) and (1,1,1) and the 3 other box diagonals (also",
300         "not independent). For each k-vector the sine and cosine are used, in",
301         "combination with the velocity in 2 perpendicular directions. This gives",
302         "a total of 16*2*2=64 transverse currents. One autocorrelation is",
303         "calculated fitted for each k-vector, which gives 16 TCAFs. Each of",
304         "these TCAFs is fitted to [MATH]f(t) = [EXP]-v[exp]([COSH]Wv[cosh] + 1/W ",
305         "[SINH]Wv[sinh])[math],",
306         "[MATH]v = -t/(2 [GRK]tau[grk])[math], [MATH]W = [SQRT]1 - 4 [GRK]tau[grk] ",
307         "[GRK]eta[grk]/[GRK]rho[grk] k^2[sqrt][math], which gives 16 values of [GRK]tau[grk]",
308         "and [GRK]eta[grk]. The fit weights decay exponentially with time constant [MATH]w[math] ",
309         "(given with [TT]-wt[tt]) as [MATH][EXP]-t/w[exp][math], and the TCAF and",
310         "fit are calculated up to time [MATH]5*w[math].",
311         "The [GRK]eta[grk] values should be fitted to [MATH]1 - a [GRK]eta[grk](k) k^2[math], ",
312         "from which one can estimate the shear viscosity at k=0.[PAR]",
313         "When the box is cubic, one can use the option [TT]-oc[tt], which",
314         "averages the TCAFs over all k-vectors with the same length.",
315         "This results in more accurate TCAFs.",
316         "Both the cubic TCAFs and fits are written to [TT]-oc[tt]",
317         "The cubic [GRK]eta[grk] estimates are also written to [TT]-ov[tt].[PAR]",
318         "With option [TT]-mol[tt], the transverse current is determined of",
319         "molecules instead of atoms. In this case, the index group should",
320         "consist of molecule numbers instead of atom numbers.[PAR]",
321         "The k-dependent viscosities in the [TT]-ov[tt] file should be",
322         "fitted to [MATH][GRK]eta[grk](k) = [GRK]eta[grk][SUB]0[sub] (1 - a k^2)[math] to obtain ",
323         "the viscosity at",
324         "infinite wavelength.[PAR]",
325         "[BB]Note:[bb] make sure you write coordinates and velocities often enough.",
326         "The initial, non-exponential, part of the autocorrelation function",
327         "is very important for obtaining a good fit."
328     };
329
330     static gmx_bool bMol = FALSE, bK34 = FALSE;
331     static real     wt   = 5;
332     t_pargs         pa[] = {
333         { "-mol", FALSE, etBOOL, { &bMol }, "Calculate TCAF of molecules" },
334         { "-k34", FALSE, etBOOL, { &bK34 }, "Also use k=(3,0,0) and k=(4,0,0)" },
335         { "-wt", FALSE, etREAL, { &wt }, "Exponential decay time for the TCAF fit weights" }
336     };
337
338     t_topology        top;
339     PbcType           pbcType;
340     t_trxframe        fr;
341     matrix            box;
342     gmx_bool          bTop;
343     int               gnx;
344     int *             index, *atndx = nullptr, at;
345     char*             grpname;
346     char              title[256];
347     real              t0, t1, dt, m, mtot, sysmass, rho, sx, cx;
348     t_trxstatus*      status;
349     int               nframes, n_alloc, i, j, k, d;
350     rvec              mv_mol, cm_mol, kfac[NK];
351     int               nkc, nk, ntc;
352     real**            tc;
353     gmx_output_env_t* oenv;
354
355 #define NHISTO 360
356
357     t_filenm fnm[] = { { efTRN, "-f", nullptr, ffREAD },      { efTPS, nullptr, nullptr, ffOPTRD },
358                        { efNDX, nullptr, nullptr, ffOPTRD },  { efXVG, "-ot", "transcur", ffOPTWR },
359                        { efXVG, "-oa", "tcaf_all", ffWRITE }, { efXVG, "-o", "tcaf", ffWRITE },
360                        { efXVG, "-of", "tcaf_fit", ffWRITE }, { efXVG, "-oc", "tcaf_cub", ffOPTWR },
361                        { efXVG, "-ov", "visc_k", ffWRITE } };
362 #define NFILE asize(fnm)
363     int      npargs;
364     t_pargs* ppa;
365
366     npargs = asize(pa);
367     ppa    = add_acf_pargs(&npargs, pa);
368
369     if (!parse_common_args(
370                 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
371     {
372         sfree(ppa);
373         return 0;
374     }
375
376     bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, nullptr, nullptr, box, TRUE);
377     get_index(&top.atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
378
379     if (bMol)
380     {
381         if (!bTop)
382         {
383             gmx_fatal(FARGS, "Need a topology to determine the molecules");
384         }
385         atndx = top.mols.index;
386     }
387
388     if (bK34)
389     {
390         nkc = NKC;
391     }
392     else
393     {
394         nkc = NKC0;
395     }
396     nk = kset_c[nkc];
397     GMX_ASSERT(nk >= 16, "Has to be over 16 because nkc is either NKC or NKC0.");
398     ntc = nk * NPK;
399
400     sprintf(title, "Velocity Autocorrelation Function for %s", grpname);
401
402     sysmass = 0;
403     for (i = 0; i < nk; i++)
404     {
405         if (iprod(v0[i], v1[i]) != 0)
406         {
407             gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
408         }
409         if (iprod(v0[i], v2[i]) != 0)
410         {
411             gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
412         }
413         if (iprod(v1[i], v2[i]) != 0)
414         {
415             gmx_fatal(FARGS, "DEATH HORROR: vectors not orthogonal");
416         }
417         unitv(v1[i], v1[i]);
418         unitv(v2[i], v2[i]);
419     }
420     snew(tc, ntc);
421     for (i = 0; i < top.atoms.nr; i++)
422     {
423         sysmass += top.atoms.atom[i].m;
424     }
425
426     read_first_frame(oenv, &status, ftp2fn(efTRN, NFILE, fnm), &fr, TRX_NEED_X | TRX_NEED_V);
427     t0 = fr.time;
428
429     n_alloc = 0;
430     nframes = 0;
431     rho     = 0;
432
433     do
434     {
435
436         if (nframes >= n_alloc)
437         {
438             n_alloc += 100;
439             for (i = 0; i < ntc; i++)
440             {
441                 srenew(tc[i], n_alloc);
442             }
443         }
444
445         rho += 1 / det(fr.box);
446         for (k = 0; k < nk; k++)
447         {
448             for (d = 0; d < DIM; d++)
449             {
450                 kfac[k][d] = 2 * M_PI * v0[k][d] / fr.box[d][d];
451             }
452         }
453         for (i = 0; i < ntc; i++)
454         {
455             tc[i][nframes] = 0;
456         }
457
458         for (i = 0; i < gnx; i++)
459         {
460             if (bMol)
461             {
462                 clear_rvec(mv_mol);
463                 clear_rvec(cm_mol);
464                 mtot = 0;
465                 for (j = 0; j < atndx[index[i] + 1] - atndx[index[i]]; j++)
466                 {
467                     at = atndx[index[i]] + j;
468                     m  = top.atoms.atom[at].m;
469                     mv_mol[XX] += m * fr.v[at][XX];
470                     mv_mol[YY] += m * fr.v[at][YY];
471                     mv_mol[ZZ] += m * fr.v[at][ZZ];
472                     cm_mol[XX] += m * fr.x[at][XX];
473                     cm_mol[YY] += m * fr.x[at][YY];
474                     cm_mol[ZZ] += m * fr.x[at][ZZ];
475                     mtot += m;
476                 }
477                 svmul(1.0 / mtot, cm_mol, cm_mol);
478             }
479             else
480             {
481                 svmul(top.atoms.atom[index[i]].m, fr.v[index[i]], mv_mol);
482             }
483
484             if (!bMol)
485             {
486                 copy_rvec(fr.x[index[i]], cm_mol);
487             }
488             j = 0;
489             for (k = 0; k < nk; k++)
490             {
491                 sx = std::sin(iprod(kfac[k], cm_mol));
492                 cx = std::cos(iprod(kfac[k], cm_mol));
493                 tc[j][nframes] += sx * iprod(v1[k], mv_mol);
494                 j++;
495                 tc[j][nframes] += cx * iprod(v1[k], mv_mol);
496                 j++;
497                 tc[j][nframes] += sx * iprod(v2[k], mv_mol);
498                 j++;
499                 tc[j][nframes] += cx * iprod(v2[k], mv_mol);
500                 j++;
501             }
502         }
503
504         t1 = fr.time;
505         nframes++;
506     } while (read_next_frame(oenv, status, &fr));
507     close_trx(status);
508
509     dt = (t1 - t0) / (nframes - 1);
510
511     rho *= sysmass / nframes * AMU / (NANO * NANO * NANO);
512     fprintf(stdout, "Density = %g (kg/m^3)\n", rho);
513     process_tcaf(nframes,
514                  dt,
515                  nkc,
516                  tc,
517                  kfac,
518                  rho,
519                  wt,
520                  opt2fn_null("-ot", NFILE, fnm),
521                  opt2fn("-oa", NFILE, fnm),
522                  opt2fn("-o", NFILE, fnm),
523                  opt2fn("-of", NFILE, fnm),
524                  opt2fn_null("-oc", NFILE, fnm),
525                  opt2fn("-ov", NFILE, fnm),
526                  oenv);
527
528     return 0;
529 }