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     "g_velacc 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."
 120   };
 121
 122   static bool bM=FALSE,bMol=FALSE;
 123   t_pargs pa[] = {
 124     { "-m", FALSE, etBOOL, {&bM},
 125       "Calculate the momentum autocorrelation function" },
 126     { "-mol", FALSE, etBOOL, {&bMol},
 127       "Calculate the velocity acf of molecules" }
 128   };
 129
 130   t_topology top;
 131   int        ePBC=-1;
 132   t_trxframe fr;
 133   matrix     box;
 134   bool       bTPS=FALSE,bTop=FALSE;
 135   int        gnx;
 136   atom_id    *index;
 137   char       *grpname;
 138   char       title[256];
 139   real       t0,t1,m;
 140   int        status,teller,n_alloc,i,j,tel3,k,l;
 141   rvec       mv_mol;
 142   real       **c1;
 143   real       *normm=NULL;
 144   output_env_t oenv;
 145
 146 #define NHISTO 360
 147
 148   t_filenm  fnm[] = {
 149     { efTRN, "-f",    NULL,   ffREAD  },
 150     { efTPS, NULL,    NULL,   ffOPTRD },
 151     { efNDX, NULL,    NULL,   ffOPTRD },
 152     { efXVG, "-o",    "vac",  ffWRITE }
 153   };
 154 #define NFILE asize(fnm)
 155   int     npargs;
 156   t_pargs *ppa;
 157
 158   CopyRight(stderr,argv[0]);
 159   npargs = asize(pa);
 160   ppa    = add_acf_pargs(&npargs,pa);
 161   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 162                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
 163
 164   if (bMol)
 165     bTPS = bM || ftp2bSet(efTPS,NFILE,fnm) || !ftp2bSet(efNDX,NFILE,fnm);
 166
 167   if (bTPS) {
 168     bTop=read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,NULL,NULL,box,
 169                        TRUE);
 170     get_index(&top.atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 171   } else
 172     rd_index(ftp2fn(efNDX,NFILE,fnm),1,&gnx,&index,&grpname);
 173
 174   if (bMol) {
 175     if (!bTop)
 176       gmx_fatal(FARGS,"Need a topology to determine the molecules");
 177     snew(normm,top.atoms.nr);
 178     precalc(top,normm);
 179     index_atom2mol(&gnx,index,&top.mols);
 180   }
 181
 182   /* Correlation stuff */
 183   snew(c1,gnx);
 184   for(i=0; (i<gnx); i++)
 185     c1[i]=NULL;
 186
 187   read_first_frame(oenv,&status,ftp2fn(efTRN,NFILE,fnm),&fr,TRX_NEED_V);
 188   t0=fr.time;
 189
 190   n_alloc=0;
 191   teller=0;
 192   do {
 193     if (teller >= n_alloc) {
 194       n_alloc+=100;
 195       for(i=0; i<gnx; i++)
 196         srenew(c1[i],DIM*n_alloc);
 197     }
 198     tel3=3*teller;
 199     if (bMol)
 200       for(i=0; i<gnx; i++) {
 201         clear_rvec(mv_mol);
 202         k=top.mols.index[index[i]];
 203         l=top.mols.index[index[i]+1];
 204         for(j=k; j<l; j++) {
 205           if (bM)
 206             m = top.atoms.atom[j].m;
 207           else
 208             m = normm[j];
 209           mv_mol[XX] += m*fr.v[j][XX];
 210           mv_mol[YY] += m*fr.v[j][YY];
 211           mv_mol[ZZ] += m*fr.v[j][ZZ];
 212         }
 213         c1[i][tel3+XX]=mv_mol[XX];
 214         c1[i][tel3+YY]=mv_mol[YY];
 215         c1[i][tel3+ZZ]=mv_mol[ZZ];
 216       }
 217      else
 218       for(i=0; i<gnx; i++) {
 219         if (bM)
 220           m = top.atoms.atom[index[i]].m;
 221         else
 222          m = 1;
 223         c1[i][tel3+XX]=m*fr.v[index[i]][XX];
 224         c1[i][tel3+YY]=m*fr.v[index[i]][YY];
 225         c1[i][tel3+ZZ]=m*fr.v[index[i]][ZZ];
 226       }
 227
 228     t1=fr.time;
 229
 230     teller ++;
 231   } while (read_next_frame(oenv,status,&fr));
 232
 233         close_trj(status);
 234
 235   do_autocorr(ftp2fn(efXVG,NFILE,fnm), oenv,
 236               bM ?
 237               "Momentum Autocorrelation Function" :
 238               "Velocity Autocorrelation Function",
 239               teller,gnx,c1,(t1-t0)/(teller-1),eacVector,TRUE);
 240
 241   do_view(oenv,ftp2fn(efXVG,NFILE,fnm),"-nxy");
 242
 243   thanx(stderr);
 244
 245   return 0;
 246 }