src/tools/gmx_tcaf.c

   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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * 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 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41 #include <stdio.h>
  42
  43 #include <math.h>
  44 #include "confio.h"
  45 #include "copyrite.h"
  46 #include "gmx_fatal.h"
  47 #include "futil.h"
  48 #include "gstat.h"
  49 #include "macros.h"
  50 #include "maths.h"
  51 #include "physics.h"
  52 #include "index.h"
  53 #include "smalloc.h"
  54 #include "statutil.h"
  55 #include <string.h>
  56 #include "sysstuff.h"
  57 #include "txtdump.h"
  58 #include "typedefs.h"
  59 #include "vec.h"
  60 #include "strdb.h"
  61 #include "xvgr.h"
  62 #include "pbc.h"
  63 #include "gmx_ana.h"
  64
  65
  66 #define NK  24
  67 #define NPK 4
  68
  69 #define NKC  6
  70 #define NKC0 4
  71 int kset_c[NKC+1] = { 0, 3, 9, 13, 16, 19, NK };
  72
  73 rvec v0[NK]={{1,0,0},{0,1,0},{0,0,1}, {1, 1,0},{1,-1,0},{1,0,1}, {1,0,-1},{0,1, 1},{0,1,-1}, {1, 1, 1},{1, 1,-1},{1,-1, 1},{-1,1, 1}, {2,0,0},{0,2,0},{0,0,2}, {3,0,0},{0,3,0},{0,0,3}, {4,0,0},{0,4,0},{0,0,4}};
  74 rvec v1[NK]={{0,1,0},{0,0,1},{1,0,0}, {0, 0,1},{0, 0,1},{0,1,0}, {0,1, 0},{1,0, 0},{1,0, 0}, {1,-1, 0},{1,-1, 0},{1, 0,-1},{ 0,1,-1}, {0,1,0},{0,0,1},{1,0,0}, {0,1,0},{0,0,1},{1,0,0}, {0,1,0},{0,0,1},{1,0,0}};
  75 rvec v2[NK]={{0,0,1},{1,0,0},{0,1,0}, {1,-1,0},{1, 1,0},{1,0,-1},{1,0, 1},{0,1,-1},{0,1, 1}, {1, 1,-2},{1, 1, 2},{1, 2, 1},{ 2,1, 1}, {0,0,1},{1,0,0},{0,1,0}, {0,0,1},{1,0,0},{0,1,0}, {0,0,1},{1,0,0},{0,1,0}};
  76
  77 static void process_tcaf(int nframes,real dt,int nkc,real **tc,rvec *kfac,
  78                          real rho,real wt,const char *fn_trans,
  79                          const char *fn_tca,const char *fn_tc,
  80                          const char *fn_tcf,const char *fn_cub,
  81                          const char *fn_vk, const output_env_t oenv)
  82 {
  83   FILE *fp,*fp_vk,*fp_cub=NULL;
  84   int  nk,ntc;
  85   real **tcaf,**tcafc=NULL,eta;
  86   int  i,j,k,kc;
  87   int  ncorr;
  88   real fitparms[3],*sig;
  89
  90   nk  = kset_c[nkc];
  91   ntc = nk*NPK;
  92
  93   if (fn_trans) {
  94     fp = xvgropen(fn_trans,"Transverse Current","Time (ps)","TC (nm/ps)",
  95                   oenv);
  96     for(i=0; i<nframes; i++) {
  97       fprintf(fp,"%g",i*dt);
  98       for(j=0; j<ntc; j++)
  99         fprintf(fp," %g",tc[j][i]);
 100       fprintf(fp,"\n");
 101     }
 102     ffclose(fp);
 103     do_view(oenv,fn_trans,"-nxy");
 104   }
 105
 106   ncorr = (nframes+1)/2;
 107   if (ncorr > (int)(5*wt/dt+0.5))
 108     ncorr = (int)(5*wt/dt+0.5)+1;
 109   snew(tcaf,nk);
 110   for(k=0; k<nk; k++)
 111     snew(tcaf[k],ncorr);
 112   if (fn_cub) {
 113      snew(tcafc,nkc);
 114      for(k=0; k<nkc; k++)
 115        snew(tcafc[k],ncorr);
 116   }
 117   snew(sig,ncorr);
 118   for(i=0; i<ncorr; i++)
 119     sig[i]=exp(0.5*i*dt/wt);
 120
 121   low_do_autocorr(fn_tca,oenv,"Transverse Current Autocorrelation Functions",
 122                   nframes,ntc,ncorr,tc,dt,eacNormal,
 123                   1,FALSE,FALSE,FALSE,0,0,0,0);
 124   do_view(oenv,fn_tca,"-nxy");
 125
 126   fp = xvgropen(fn_tc,"Transverse Current Autocorrelation Functions",
 127                 "Time (ps)","TCAF",oenv);
 128   for(i=0; i<ncorr; i++) {
 129     kc = 0;
 130     fprintf(fp,"%g",i*dt);
 131     for(k=0; k<nk; k++) {
 132       for(j=0; j<NPK; j++)
 133         tcaf[k][i] += tc[NPK*k+j][i];
 134       if (fn_cub)
 135         for(j=0; j<NPK; j++)
 136           tcafc[kc][i] += tc[NPK*k+j][i];
 137       if (i == 0)
 138         fprintf(fp," %g",1.0);
 139       else {
 140         tcaf[k][i] /= tcaf[k][0];
 141         fprintf(fp," %g",tcaf[k][i]);
 142       }
 143       if (k+1 == kset_c[kc+1])
 144         kc++;
 145     }
 146     fprintf(fp,"\n");
 147   }
 148   ffclose(fp);
 149   do_view(oenv,fn_tc,"-nxy");
 150
 151   if (fn_cub) {
 152     fp_cub = xvgropen(fn_cub,"TCAFs and fits", "Time (ps)","TCAF",oenv);
 153     for(kc=0; kc<nkc; kc++) {
 154       fprintf(fp_cub,"%g %g\n",0.0,1.0);
 155       for(i=1; i<ncorr; i++) {
 156         tcafc[kc][i] /= tcafc[kc][0];
 157         fprintf(fp_cub,"%g %g\n",i*dt,tcafc[kc][i]);
 158       }
 159       fprintf(fp_cub,"&\n");
 160       tcafc[kc][0] = 1.0;
 161     }
 162   }
 163
 164   fp_vk = xvgropen(fn_vk,"Fits","k (nm\\S-1\\N)",
 165                    "\\8h\\4 (10\\S-3\\N kg m\\S-1\\N s\\S-1\\N)",oenv);
 166   fprintf(fp_vk,"@    s0 symbol 2\n");
 167   fprintf(fp_vk,"@    s0 symbol color 1\n");
 168   fprintf(fp_vk,"@    s0 linestyle 0\n");
 169   if (fn_cub) {
 170     fprintf(fp_vk,"@    s1 symbol 3\n");
 171     fprintf(fp_vk,"@    s1 symbol color 2\n");
 172   }
 173   fp = xvgropen(fn_tcf,"TCAF Fits","Time (ps)","",oenv);
 174   for(k=0; k<nk; k++) {
 175     tcaf[k][0] = 1.0;
 176     fitparms[0]  = 1;
 177     fitparms[1]  = 1;
 178     do_lmfit(ncorr,tcaf[k],sig,dt,0,0,ncorr*dt,
 179              oenv,bDebugMode(),effnVAC,fitparms,0);
 180     eta = 1000*fitparms[1]*rho/
 181       (4*fitparms[0]*PICO*norm2(kfac[k])/(NANO*NANO));
 182     fprintf(stdout,"k %6.3f  tau %6.3f  eta %8.5f 10^-3 kg/(m s)\n",
 183             norm(kfac[k]),fitparms[0],eta);
 184     fprintf(fp_vk,"%6.3f %g\n",norm(kfac[k]),eta);
 185     for(i=0; i<ncorr; i++)
 186       fprintf(fp,"%g %g\n",i*dt,fit_function(effnVAC,fitparms,i*dt));
 187     fprintf(fp,"&\n");
 188   }
 189   ffclose(fp);
 190   do_view(oenv,fn_tcf,"-nxy");
 191
 192   if (fn_cub) {
 193     fprintf(stdout,"Averaged over k-vectors:\n");
 194     fprintf(fp_vk,"&\n");
 195     for(k=0; k<nkc; k++) {
 196       tcafc[k][0] = 1.0;
 197       fitparms[0]  = 1;
 198       fitparms[1]  = 1;
 199       do_lmfit(ncorr,tcafc[k],sig,dt,0,0,ncorr*dt,
 200                oenv,bDebugMode(),effnVAC,fitparms,0);
 201       eta = 1000*fitparms[1]*rho/
 202         (4*fitparms[0]*PICO*norm2(kfac[kset_c[k]])/(NANO*NANO));
 203       fprintf(stdout,
 204               "k %6.3f  tau %6.3f  Omega %6.3f  eta %8.5f 10^-3 kg/(m s)\n",
 205               norm(kfac[kset_c[k]]),fitparms[0],fitparms[1],eta);
 206       fprintf(fp_vk,"%6.3f %g\n",norm(kfac[kset_c[k]]),eta);
 207       for(i=0; i<ncorr; i++)
 208         fprintf(fp_cub,"%g %g\n",i*dt,fit_function(effnVAC,fitparms,i*dt));
 209       fprintf(fp_cub,"&\n");
 210     }
 211     fprintf(fp_vk,"&\n");
 212     ffclose(fp_cub);
 213     do_view(oenv,fn_cub,"-nxy");
 214   }
 215   ffclose(fp_vk);
 216   do_view(oenv,fn_vk,"-nxy");
 217 }
 218
 219
 220 int gmx_tcaf(int argc,char *argv[])
 221 {
 222   const char *desc[] = {
 223     "[TT]g_tcaf[tt] computes tranverse current autocorrelations.",
 224     "These are used to estimate the shear viscosity, [GRK]eta[grk].",
 225     "For details see: Palmer, Phys. Rev. E 49 (1994) pp 359-366.[PAR]",
 226     "Transverse currents are calculated using the",
 227     "k-vectors (1,0,0) and (2,0,0) each also in the [IT]y[it]- and [IT]z[it]-direction,",
 228     "(1,1,0) and (1,-1,0) each also in the 2 other planes (these vectors",
 229     "are not independent) and (1,1,1) and the 3 other box diagonals (also",
 230     "not independent). For each k-vector the sine and cosine are used, in",
 231     "combination with the velocity in 2 perpendicular directions. This gives",
 232     "a total of 16*2*2=64 transverse currents. One autocorrelation is",
 233     "calculated fitted for each k-vector, which gives 16 TCAFs. Each of",
 234     "these TCAFs is fitted to [MATH]f(t) = [EXP]-v[exp]([COSH]Wv[cosh] + 1/W [SINH]Wv[sinh])[math],",
 235     "[MATH]v = -t/(2 [GRK]tau[grk])[math], [MATH]W = [SQRT]1 - 4 [GRK]tau[grk] [GRK]eta[grk]/[GRK]rho[grk] k^2[sqrt][math], which gives 16 values of [GRK]tau[grk]",
 236     "and [GRK]eta[grk]. The fit weights decay exponentially with time constant [MATH]w[math] (given with [TT]-wt[tt]) as [MATH][EXP]-t/w[exp][math], and the TCAF and",
 237     "fit are calculated up to time [MATH]5*w[math].",
 238     "The [GRK]eta[grk] values should be fitted to [MATH]1 - a [GRK]eta[grk](k) k^2[math], from which",
 239     "one can estimate the shear viscosity at k=0.[PAR]",
 240     "When the box is cubic, one can use the option [TT]-oc[tt], which",
 241     "averages the TCAFs over all k-vectors with the same length.",
 242     "This results in more accurate TCAFs.",
 243     "Both the cubic TCAFs and fits are written to [TT]-oc[tt]",
 244     "The cubic [GRK]eta[grk] estimates are also written to [TT]-ov[tt].[PAR]",
 245     "With option [TT]-mol[tt], the transverse current is determined of",
 246     "molecules instead of atoms. In this case, the index group should",
 247     "consist of molecule numbers instead of atom numbers.[PAR]",
 248     "The k-dependent viscosities in the [TT]-ov[tt] file should be",
 249     "fitted to [MATH][GRK]eta[grk](k) = [GRK]eta[grk][SUB]0[sub] (1 - a k^2)[math] to obtain the viscosity at",
 250     "infinite wavelength.[PAR]",
 251     "[BB]Note:[bb] make sure you write coordinates and velocities often enough.",
 252     "The initial, non-exponential, part of the autocorrelation function",
 253     "is very important for obtaining a good fit."
 254   };
 255
 256   static gmx_bool bMol=FALSE,bK34=FALSE;
 257   static real wt=5;
 258   t_pargs pa[] = {
 259     { "-mol", FALSE, etBOOL, {&bMol},
 260       "Calculate TCAF of molecules" },
 261     { "-k34", FALSE, etBOOL, {&bK34},
 262       "Also use k=(3,0,0) and k=(4,0,0)" },
 263     { "-wt", FALSE, etREAL, {&wt},
 264       "Exponential decay time for the TCAF fit weights" }
 265   };
 266
 267   t_topology top;
 268   int        ePBC;
 269   t_trxframe fr;
 270   matrix     box;
 271   gmx_bool       bTPS,bTop; /* ,bCubic; */
 272   int        gnx;
 273   atom_id    *index,*atndx=NULL,at;
 274   char       *grpname;
 275   char       title[256];
 276   real       t0,t1,dt,m,mtot,sysmass,rho,sx,cx;
 277   t_trxstatus *status;
 278   int        nframes,n_alloc,i,j,k,d;
 279   rvec       mv_mol,cm_mol,kfac[NK];
 280   int        nkc,nk,ntc;
 281   real       **c1,**tc;
 282   output_env_t oenv;
 283
 284 #define NHISTO 360
 285
 286   t_filenm  fnm[] = {
 287     { efTRN, "-f",    NULL,      ffREAD  },
 288     { efTPS, NULL,    NULL,      ffOPTRD },
 289     { efNDX, NULL,    NULL,      ffOPTRD },
 290     { efXVG, "-ot",  "transcur", ffOPTWR },
 291     { efXVG, "-oa",  "tcaf_all", ffWRITE },
 292     { efXVG, "-o",   "tcaf",     ffWRITE },
 293     { efXVG, "-of",  "tcaf_fit", ffWRITE },
 294     { efXVG, "-oc",  "tcaf_cub", ffOPTWR },
 295     { efXVG, "-ov",  "visc_k",   ffWRITE }
 296   };
 297 #define NFILE asize(fnm)
 298   int     npargs;
 299   t_pargs *ppa;
 300
 301   CopyRight(stderr,argv[0]);
 302   npargs = asize(pa);
 303   ppa    = add_acf_pargs(&npargs,pa);
 304
 305   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 306                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
 307
 308   bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,NULL,NULL,box,
 309                      TRUE);
 310   get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 311
 312   if (bMol) {
 313     if (!bTop)
 314       gmx_fatal(FARGS,"Need a topology to determine the molecules");
 315     atndx = top.mols.index;
 316   }
 317
 318   if (bK34)
 319     nkc = NKC;
 320   else
 321     nkc = NKC0;
 322   nk  = kset_c[nkc];
 323   ntc = nk*NPK;
 324
 325   sprintf(title,"Velocity Autocorrelation Function for %s",grpname);
 326
 327   sysmass = 0;
 328   for(i=0; i<nk; i++) {
 329     if (iprod(v0[i],v1[i]) != 0)
 330       gmx_fatal(FARGS,"DEATH HORROR: vectors not orthogonal");
 331     if (iprod(v0[i],v2[i]) != 0)
 332       gmx_fatal(FARGS,"DEATH HORROR: vectors not orthogonal");
 333     if (iprod(v1[i],v2[i]) != 0)
 334         gmx_fatal(FARGS,"DEATH HORROR: vectors not orthogonal");
 335     unitv(v1[i],v1[i]);
 336     unitv(v2[i],v2[i]);
 337   }
 338   snew(tc,ntc);
 339     for(i=0; i<top.atoms.nr; i++)
 340       sysmass += top.atoms.atom[i].m;
 341
 342   read_first_frame(oenv,&status,ftp2fn(efTRN,NFILE,fnm),&fr,
 343                    TRX_NEED_X | TRX_NEED_V);
 344   t0=fr.time;
 345
 346   n_alloc=0;
 347   nframes=0;
 348   rho=0;
 349   /* bCubic = TRUE; */
 350   do {
 351     /*
 352     bCubic = bCubic && !TRICLINIC(fr.box) &&
 353     fabs(fr.box[XX][XX]-fr.box[YY][YY]) < 0.001*fr.box[XX][XX] &&
 354       fabs(fr.box[XX][XX]-fr.box[ZZ][ZZ]) < 0.001*fr.box[XX][XX];
 355       */
 356
 357     if (nframes >= n_alloc) {
 358       n_alloc+=100;
 359       for (i=0; i<ntc; i++)
 360         srenew(tc[i],n_alloc);
 361     }
 362
 363     rho += 1/det(fr.box);
 364     for(k=0; k<nk; k++)
 365       for(d=0; d<DIM; d++)
 366         kfac[k][d] = 2*M_PI*v0[k][d]/fr.box[d][d];
 367     for (i=0; i<ntc; i++)
 368       tc[i][nframes] = 0;
 369
 370     for(i=0; i<gnx; i++) {
 371       if (bMol) {
 372         clear_rvec(mv_mol);
 373         clear_rvec(cm_mol);
 374         mtot = 0;
 375         for(j=0; j<atndx[index[i]+1] - atndx[index[i]]; j++) {
 376           at = atndx[index[i]] + j;
 377           m  = top.atoms.atom[at].m;
 378           mv_mol[XX] += m*fr.v[at][XX];
 379           mv_mol[YY] += m*fr.v[at][YY];
 380           mv_mol[ZZ] += m*fr.v[at][ZZ];
 381           cm_mol[XX] += m*fr.x[at][XX];
 382           cm_mol[YY] += m*fr.x[at][YY];
 383           cm_mol[ZZ] += m*fr.x[at][ZZ];
 384           mtot += m;
 385         }
 386         svmul(1.0/mtot,cm_mol,cm_mol);
 387       } else
 388         svmul(top.atoms.atom[index[i]].m,fr.v[index[i]],mv_mol);
 389
 390       if (!bMol)
 391         copy_rvec(fr.x[index[i]],cm_mol);
 392       j=0;
 393       for(k=0; k<nk; k++) {
 394         sx = sin(iprod(kfac[k],cm_mol));
 395         cx = cos(iprod(kfac[k],cm_mol));
 396         tc[j][nframes] += sx*iprod(v1[k],mv_mol);
 397         j++;
 398         tc[j][nframes] += cx*iprod(v1[k],mv_mol);
 399         j++;
 400         tc[j][nframes] += sx*iprod(v2[k],mv_mol);
 401         j++;
 402         tc[j][nframes] += cx*iprod(v2[k],mv_mol);
 403         j++;
 404       }
 405     }
 406
 407     t1=fr.time;
 408     nframes ++;
 409   } while (read_next_frame(oenv,status,&fr));
 410   close_trj(status);
 411
 412   dt = (t1-t0)/(nframes-1);
 413
 414   rho *= sysmass/nframes*AMU/(NANO*NANO*NANO);
 415   fprintf(stdout,"Density = %g (kg/m^3)\n",rho);
 416   process_tcaf(nframes,dt,nkc,tc,kfac,rho,wt,
 417                opt2fn_null("-ot",NFILE,fnm),
 418                opt2fn("-oa",NFILE,fnm),opt2fn("-o",NFILE,fnm),
 419                opt2fn("-of",NFILE,fnm),opt2fn_null("-oc",NFILE,fnm),
 420                opt2fn("-ov",NFILE,fnm),oenv);
 421
 422   thanx(stderr);
 423
 424   return 0;
 425 }