src/tools/gmx_velacc.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 #include <math.h>
  40
  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 "gmx_ana.h"
  60
  61
  62 static void index_atom2mol(int *n,atom_id *index,t_block *mols)
  63 {
  64   int nat,i,nmol,mol,j;
  65
  66   nat = *n;
  67   i = 0;
  68   nmol = 0;
  69   mol = 0;
  70   while (i < nat) {
  71     while (index[i] > mols->index[mol]) {
  72       mol++;
  73       if (mol >= mols->nr)
  74         gmx_fatal(FARGS,"Atom index out of range: %d",index[i]+1);
  75     }
  76     for(j=mols->index[mol]; j<mols->index[mol+1]; j++) {
  77       if (i >= nat || index[i] != j)
  78         gmx_fatal(FARGS,"The index group does not consist of whole molecules");
  79       i++;
  80     }
  81     index[nmol++] = mol;
  82   }
  83
  84   fprintf(stderr,"\nSplit group of %d atoms into %d molecules\n",nat,nmol);
  85
  86   *n = nmol;
  87 }
  88
  89 static void precalc(t_topology top,real normm[]){
  90
  91   real mtot;
  92   int i,j,k,l;
  93
  94   for(i=0;i<top.mols.nr;i++){
  95     k=top.mols.index[i];
  96     l=top.mols.index[i+1];
  97     mtot=0.0;
  98
  99     for(j=k;j<l;j++)
 100       mtot+=top.atoms.atom[j].m;
 101
 102     for(j=k;j<l;j++)
 103       normm[j]=top.atoms.atom[j].m/mtot;
 104
 105   }
 106
 107 }
 108
 109
 110
 111 int gmx_velacc(int argc,char *argv[])
 112 {
 113   const char *desc[] = {
 114     "[TT]g_velacc[tt] computes the velocity autocorrelation function.",
 115     "When the [TT]-m[tt] option is used, the momentum autocorrelation",
 116     "function is calculated.[PAR]",
 117     "With option [TT]-mol[tt] the velocity autocorrelation function of",
 118     "molecules is calculated. In this case the index group should consist",
 119     "of molecule numbers instead of atom numbers.[PAR]",
 120     "Be sure that your trajectory contains frames with velocity information",
 121     "(i.e. [TT]nstvout[tt] was set in your original [TT].mdp[tt] file),",
 122     "and that the time interval between data collection points is",
 123     "much shorter than the time scale of the autocorrelation."
 124   };
 125
 126   static gmx_bool bM=FALSE,bMol=FALSE;
 127   t_pargs pa[] = {
 128     { "-m", FALSE, etBOOL, {&bM},
 129       "Calculate the momentum autocorrelation function" },
 130     { "-mol", FALSE, etBOOL, {&bMol},
 131       "Calculate the velocity acf of molecules" }
 132   };
 133
 134   t_topology top;
 135   int        ePBC=-1;
 136   t_trxframe fr;
 137   matrix     box;
 138   gmx_bool       bTPS=FALSE,bTop=FALSE;
 139   int        gnx;
 140   atom_id    *index;
 141   char       *grpname;
 142   char       title[256];
 143   real       t0,t1,m;
 144   t_trxstatus *status;
 145   int        teller,n_alloc,i,j,tel3,k,l;
 146   rvec       mv_mol;
 147   real       **c1;
 148   real       *normm=NULL;
 149   output_env_t oenv;
 150
 151 #define NHISTO 360
 152
 153   t_filenm  fnm[] = {
 154     { efTRN, "-f",    NULL,   ffREAD  },
 155     { efTPS, NULL,    NULL,   ffOPTRD },
 156     { efNDX, NULL,    NULL,   ffOPTRD },
 157     { efXVG, "-o",    "vac",  ffWRITE }
 158   };
 159 #define NFILE asize(fnm)
 160   int     npargs;
 161   t_pargs *ppa;
 162
 163   CopyRight(stderr,argv[0]);
 164   npargs = asize(pa);
 165   ppa    = add_acf_pargs(&npargs,pa);
 166   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 167                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
 168
 169   if (bMol || bM) {
 170     bTPS = ftp2bSet(efTPS,NFILE,fnm) || !ftp2bSet(efNDX,NFILE,fnm);
 171   }
 172
 173   if (bTPS) {
 174     bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,NULL,NULL,box,
 175                        TRUE);
 176     get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 177   } else
 178     rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 179
 180   if (bMol) {
 181     if (!bTop)
 182       gmx_fatal(FARGS,"Need a topology to determine the molecules");
 183     snew(normm,top.atoms.nr);
 184     precalc(top,normm);
 185     index_atom2mol(&gnx,index,&top.mols);
 186   }
 187
 188   /* Correlation stuff */
 189   snew(c1,gnx);
 190   for(i=0; (i<gnx); i++)
 191     c1[i]=NULL;
 192
 193   read_first_frame(oenv,&status,ftp2fn(efTRN,NFILE,fnm),&fr,TRX_NEED_V);
 194   t0=fr.time;
 195
 196   n_alloc=0;
 197   teller=0;
 198   do {
 199     if (teller >= n_alloc) {
 200       n_alloc+=100;
 201       for(i=0; i<gnx; i++)
 202         srenew(c1[i],DIM*n_alloc);
 203     }
 204     tel3=3*teller;
 205     if (bMol)
 206       for(i=0; i<gnx; i++) {
 207         clear_rvec(mv_mol);
 208         k=top.mols.index[index[i]];
 209         l=top.mols.index[index[i]+1];
 210         for(j=k; j<l; j++) {
 211           if (bM)
 212             m = top.atoms.atom[j].m;
 213           else
 214             m = normm[j];
 215           mv_mol[XX] += m*fr.v[j][XX];
 216           mv_mol[YY] += m*fr.v[j][YY];
 217           mv_mol[ZZ] += m*fr.v[j][ZZ];
 218         }
 219         c1[i][tel3+XX]=mv_mol[XX];
 220         c1[i][tel3+YY]=mv_mol[YY];
 221         c1[i][tel3+ZZ]=mv_mol[ZZ];
 222       }
 223      else
 224       for(i=0; i<gnx; i++) {
 225         if (bM)
 226           m = top.atoms.atom[index[i]].m;
 227         else
 228          m = 1;
 229         c1[i][tel3+XX]=m*fr.v[index[i]][XX];
 230         c1[i][tel3+YY]=m*fr.v[index[i]][YY];
 231         c1[i][tel3+ZZ]=m*fr.v[index[i]][ZZ];
 232       }
 233
 234     t1=fr.time;
 235
 236     teller ++;
 237   } while (read_next_frame(oenv,status,&fr));
 238
 239         close_trj(status);
 240
 241   do_autocorr(ftp2fn(efXVG,NFILE,fnm), oenv,
 242               bM ?
 243               "Momentum Autocorrelation Function" :
 244               "Velocity Autocorrelation Function",
 245               teller,gnx,c1,(t1-t0)/(teller-1),eacVector,TRUE);
 246
 247   do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");
 248
 249   thanx(stderr);
 250
 251   return 0;
 252 }