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