src/tools/g_diffuse.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_diffuse_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <string.h>
  33 #include <ctype.h>
  34 #include <math.h>
  35 #include "sysstuff.h"
  36 #include "smalloc.h"
  37 #include "macros.h"
  38 #include "statutil.h"
  39 #include "maths.h"
  40 #include "futil.h"
  41 #include "vec.h"
  42 #include "rdgroup.h"
  43 #include "copyrite.h"
  44 #include "typedefs.h"
  45 #include "xvgr.h"
  46 #include "gstat.h"
  47
  48 void prepare_data(int gnx,atom_id index[],rvec xcur[],rvec xprev[],
  49                   matrix box)
  50 {
  51   int  i,m,ind;
  52   rvec hbox;
  53
  54   /* Remove periodicity */
  55   for(m=0; (m<DIM); m++)
  56     hbox[m]=0.5*box[m][m];
  57   for(i=0; (i<gnx); i++) {
  58     ind=index[i];
  59     for(m=0; (m<DIM); m++) {
  60       while(xcur[ind][m]-xprev[ind][m] <= hbox[m])
  61         xcur[ind][m] += box[m][m];
  62       while(xcur[ind][m]-xprev[ind][m] >  hbox[m])
  63         xcur[ind][m] -= box[m][m];
  64     }
  65   }
  66 }
  67
  68 static real calc_msd(rvec xref[],int nx,atom_id index[],rvec x[])
  69 {
  70   int  i,m,ind;
  71   rvec dx;
  72   real g;
  73
  74   g = 0;
  75   for(i=0; (i<nx); i++) {
  76     ind = index[i];
  77     rvec_sub(xref[ind],x[ind],dx);
  78     g += iprod(dx,dx);
  79   }
  80   return g/nx;
  81 }
  82
  83 static void analyse_msd(char *xvg,char *grpname,
  84                         int nframes,int nstart,real **msd,real xtime[],
  85                         real bfit,real efit)
  86 {
  87   FILE  *fp;
  88   char  buf[256];
  89   int   i,j,nj,minj,legind;
  90   real  *D,*Dsq,t0,msd0,D2,DeltaD;
  91   char  **leg;
  92
  93   snew(D,nstart+1);
  94   snew(Dsq,nstart+1);
  95   snew(leg,nstart+1);
  96
  97   sprintf(buf,"Diffusion constant for %s",grpname);
  98   fp=xvgropen(xvg,buf,"Time (ps)","D (cm\\S2\\Ns\\S-1\\N)");
  99   minj = nframes;
 100   for(i=0; (i<nstart); i++) {
 101     D[i]   = 0;
 102     Dsq[i] = 0;
 103     nj     = 0;
 104     for(j=0; (j<nframes); j++) {
 105       if (j > 0) {
 106         msd[i][j] *= 1000.0/(6*(xtime[j] - xtime[0]));
 107         if (msd[i][j] == 0.0)
 108           break;
 109       }
 110       if ((xtime[j] >= bfit) && (xtime[j] <= efit)) {
 111         D[i]   += msd[i][j];
 112         Dsq[i] += sqr(msd[i][j]);
 113         nj++;
 114       }
 115     }
 116     legind = (nstart > 1) ? i+1 : 0;
 117     if (nj > 0) {
 118       D[i] /= (real) nj;
 119       sprintf(buf,"D%d = %.2f, N=%d",i+1,D[i],nj);
 120       leg[legind] = strdup(buf);
 121     }
 122     else {
 123       leg[legind] = strdup("No data!");
 124     }
 125     minj=min(minj,j);
 126   }
 127   /* Compute average D */
 128   if (nstart > 1) {
 129     D[nstart] = 0.0;
 130     D2        = 0.0;
 131     for(i=0; (i<nstart); i++) {
 132       D[nstart] += D[i];
 133       D2 += sqr(D[i]);
 134     }
 135     D[nstart] /= nstart;
 136     D2        /= nstart;
 137     DeltaD = sqrt(D2 - sqr(D[nstart]));
 138     fprintf(stderr,"Dav = %.2f, RMS = %.2f\n",D[nstart],DeltaD);
 139     fprintf(fp,"# Dav = %.2f, RMS = %.2f nstart = %d\n",D[nstart],DeltaD,nstart);
 140
 141     sprintf(buf,"D\\sav\\N = %.2f",D[nstart]);
 142     leg[0] = strdup(buf);
 143   }
 144   xvgr_legend(fp,((nstart > 1) ? nstart+1 : 1),leg);
 145
 146   for(j=0; (j<minj); j++) {
 147     fprintf(fp,"%10g",xtime[j]-xtime[0]);
 148     if (nstart > 1) {
 149       msd0 = 0;
 150       for(i=0; (i<nstart); i++)
 151         msd0 += msd[i][j];
 152       fprintf(fp,"  %10g",msd0/nstart);
 153     }
 154     for(i=0; (i<nstart); i++)
 155       fprintf(fp,"  %10g",msd[i][j]);
 156     fprintf(fp,"\n");
 157   }
 158   fclose(fp);
 159 }
 160
 161 static void do_msd(char *trx,char *xvg,
 162                    int nframes,int nstart,real dt,
 163                    int nx,atom_id index[],char *grpname,
 164                    real bfit,real efit)
 165 {
 166   real   **msd,*xtime,t,t0,msdt;
 167   int    *n_offs;
 168   int    i,n,status,iframe,istart,natoms;
 169   bool   bEOF;
 170   rvec   *x[2],**xref=NULL;
 171   int    cur=0;
 172 #define prev 1-cur
 173   matrix box;
 174
 175   if (nstart <= 0)
 176     fatal_error(0,"nstart should be > 0 (iso %d)",nstart);
 177   if (nframes <= 0)
 178     fatal_error(0,"nframes should be > 0 (iso %d)",nframes);
 179
 180   natoms=read_first_x(&status,trx,&t0,&(x[cur]),box);
 181   t = t0;
 182   snew(x[prev],natoms);
 183   memcpy(x[prev],x[cur],natoms*sizeof(x[cur][0]));
 184
 185   snew(n_offs,nstart);
 186   snew(msd,nstart);
 187   snew(xref,nstart);
 188   for(i=0; (i<nstart); i++) {
 189     n_offs[i] = -1;
 190     snew(msd[i],nframes);
 191     snew(xref[i],natoms);
 192   }
 193   snew(xtime,nframes);
 194
 195   /* Check index */
 196   for(i=0; (i<nx); i++)
 197     if (index[i] >= natoms)
 198       fatal_error(0,"Index[%d] = %d should be >=0 and < %d",i,index[i],natoms);
 199
 200   /* Index is OK if we got here */
 201   iframe = 0;
 202   istart = 0;
 203   do {
 204     xtime[iframe] = t;
 205     /* Check for new starting point */
 206     if (istart < nstart) {
 207       if ((t >= (t0+istart*dt)) && (n_offs[istart] == -1)) {
 208         fprintf(stderr,"  New starting point\n");
 209         memcpy(xref[istart],x[cur],natoms*sizeof(x[cur][0]));
 210         n_offs[istart]=iframe;
 211         istart++;
 212       }
 213     }
 214     for(n=0; (n<istart); n++) {
 215       /* Compute the MSD for all starting points */
 216       msdt = calc_msd(xref[n],nx,index,x[cur]);
 217       msd[n][iframe-n_offs[n]] = msdt;
 218     }
 219
 220     /* Read new data into new array and remove PBC */
 221     cur  = prev;
 222     bEOF = !read_next_x(status,&t,natoms,x[cur],box);
 223     prepare_data(nx,index,x[cur],x[prev],box);
 224
 225     iframe++;
 226   } while ((iframe < nframes) && !bEOF);
 227   close_trj(status);
 228
 229   analyse_msd(xvg,grpname,iframe,istart,msd,xtime,bfit,efit);
 230
 231   for(i=0; (i<nstart); i++) {
 232     sfree(xref[i]);
 233     sfree(msd[i]);
 234   }
 235   sfree(xref);
 236   sfree(msd);
 237   sfree(xtime);
 238   sfree(n_offs);
 239 }
 240
 241 int main(int argc,char *argv[])
 242 {
 243   static char *desc[] = {
 244     "g_diffuse computes the diffusion constant using the mean square",
 245     "displacement of atoms from",
 246     "their initial positions (using the Einstein relation).[PAR]",
 247   };
 248   static char *bugs[] = {
 249     "When the number of starting points is too high, average and standard deviation may be meaningless"
 250   };
 251   static int  nrframes = 10;
 252   static int  nstart   = 1;
 253   static real dt       = 10.0;
 254   static real bfit=0,efit=10000;
 255   t_pargs pa[] = {
 256     { "-nframes", FALSE, etINT, &nrframes,
 257       "Number of frames in your trajectory" },
 258     { "-nstart",  FALSE, etINT, &nstart,
 259       "Number of starting points" },
 260     { "-dt",      FALSE, etREAL, &dt,
 261       "Time between starting points" },
 262     { "-beginfit",FALSE, etREAL, &bfit,
 263       "Begin fitting to a straight line at this time (ps)" },
 264     { "-endfit",  FALSE, etREAL, &efit,
 265       "End fitting here" }
 266   };
 267   t_filenm fnm[] = {
 268     { efTRX, "-f", NULL,  ffREAD },
 269     { efNDX, NULL, NULL,  ffREAD },
 270     { efXVG, NULL, "diff.xvg",  ffWRITE },
 271   };
 272 #define NFILE asize(fnm)
 273   int     nx;
 274   atom_id *index=NULL;
 275   char    *grpname=NULL;
 276
 277   CopyRight(stdout,argv[0]);
 278
 279   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 280                     NFILE,fnm,asize(pa),pa,asize(desc),desc,asize(bugs),bugs);
 281
 282   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&nx,&index,&grpname);
 283
 284   do_msd(ftp2fn(efTRX,NFILE,fnm),ftp2fn(efXVG,NFILE,fnm),
 285          nrframes,nstart,dt,nx,index,grpname,bfit,efit);
 286
 287   if (bDoView())
 288     xvgr_file(ftp2fn(efXVG,NFILE,fnm),"-nxy");
 289
 290   thanx(stdout);
 291
 292   return 0;
 293 }