src/tools/gmx_tcaf.c

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