src/tools/g_run_rms.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  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_g_run_rms_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "smalloc.h"
  33 #include "typedefs.h"
  34 #include "copyrite.h"
  35 #include "statutil.h"
  36 #include "string2.h"
  37 #include "vec.h"
  38 #include "rdgroup.h"
  39 #include "pbc.h"
  40 #include "macros.h"
  41 #include "gstat.h"
  42 #include "xvgr.h"
  43 #include "tpxio.h"
  44
  45 real calc_ener(int natoms,real w_rms[],rvec x[],rvec xp[])
  46 {
  47   real e,tmas;
  48   int  j,m;
  49
  50   /*calculate energy */
  51   e=0;
  52   tmas=0;
  53   for(j=0;(j<natoms);j++) {
  54     tmas+=w_rms[j];
  55     for(m=0;(m<3);m++)
  56       e+=w_rms[j]*(x[j][m]-xp[j][m])*(x[j][m]-xp[j][m]);
  57   }
  58   /*return energy*/
  59   return(sqrt(e/tmas));
  60 }
  61
  62 int main (int argc,char *argv[])
  63 {
  64   static char *desc[] = {
  65     "g_run_rms computes the root mean square deviation (RMSD) of a structure",
  66     "from a trajectory x(t), with respect to a reference structure from the",
  67     "same trajectory, but at a specified time before the current time",
  68     "x(t-dt) by LSQ fitting the structures on top of each other[PAR]",
  69     "This tells you something about the mobility of structure as a function",
  70     "of time"
  71   };
  72   static int run_time = 5;
  73   t_pargs pa[] = {
  74     { "-t", FALSE, etINT, &run_time,
  75       "Time interval dt between reference and actual structure" }
  76   };
  77   int          step,nre,natom,i,j,m,teller=0;
  78   real         t,lambda,*w_rls,*w_rms,tmas;
  79   int          status;
  80   t_tpxheader  header;
  81   t_inputrec   ir;
  82   t_topology   top;
  83   matrix       box;
  84   rvec         **x,*xp,*v,xcm,*temp;
  85   FILE         *fp;
  86
  87 #define RLS 0
  88 #define RMS 1
  89   int          isize[2];
  90   atom_id      *index[2];
  91   char         *grpnames[2];
  92   t_filenm     fnm[] = {
  93     { efTRX, "-f", NULL, ffREAD },
  94     { efTPX, NULL, NULL, ffREAD },
  95     { efNDX, NULL, NULL, ffREAD },
  96     { efXVG, NULL, "runrms", ffWRITE }
  97   };
  98 #define NFILE asize(fnm)
  99
 100   CopyRight(stderr,argv[0]);
 101   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 102                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 103   snew(x,run_time+1);
 104
 105   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&header);
 106   snew(xp,header.natoms);
 107   for(i=0;(i<run_time+1);i++)
 108     snew(x[i],header.natoms);
 109   snew(w_rls,header.natoms);
 110   snew(w_rms,header.natoms);
 111   snew(v,header.natoms);
 112
 113   read_tpx(ftp2fn(efTPX,NFILE,fnm),
 114            &step,&t,&lambda,&ir,
 115            box,&natom,xp,NULL,NULL,&top);
 116
 117   /*set box type*/
 118   init_pbc(box,FALSE);
 119
 120   fprintf(stderr,"select group for root least squares fit\n");
 121   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize[RLS],&index[RLS],&grpnames[RLS]);
 122   for(i=0;(i<isize[RLS]);i++)
 123     w_rls[index[RLS][i]]=top.atoms.atom[index[RLS][i]].m;
 124
 125   fprintf(stderr,"select group for root mean square calculation\n");
 126   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize[RMS],&index[RMS],&grpnames[RMS]);
 127   for(i=0;(i<isize[RMS]);i++)
 128     w_rms[index[RMS][i]]=top.atoms.atom[index[RMS][i]].m;
 129
 130   fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),"Running RMSD","Time (ps)","RMSD (nm)");
 131
 132   rm_pbc(&(top.idef),top.atoms.nr,box,xp,xp);
 133
 134   /*calculate mass*/
 135   tmas=0;
 136   for(i=0;(i<header.natoms);i++)
 137     tmas+=w_rls[i];
 138
 139   /*set center of mass to zero for reference coordinates*/
 140
 141   /*do a first step*/
 142   natom=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x[0],box);
 143   for(i=0; (i<run_time); i++) {
 144     read_next_x(status,&t,natom,x[i],box);
 145     rm_pbc(&(top.idef),top.atoms.nr,box,x[i],x[i]);
 146     for(m=0;(m<3);m++)
 147       xcm[m]=0;
 148     for(j=0;(j<header.natoms);j++)
 149       for(m=0;(m<3);m++) {
 150         xcm[m]+=x[i][j][m]*w_rls[j];
 151       }
 152     for(m=0;(m<3);m++)
 153       xcm[m]/=tmas;
 154     for(j=0;(j<header.natoms);j++) {
 155       for(m=0;(m<3);m++)
 156         x[i][j][m]-=xcm[m];
 157     }
 158   }
 159
 160   read_next_x(status,&t,natom,x[run_time],box);
 161   do {
 162     teller++;
 163     rm_pbc(&(top.idef),top.atoms.nr,box,x[run_time],x[run_time]);
 164     for(m=0;(m<3);m++)
 165       xcm[m]=0;
 166     for(j=0;(j<header.natoms);j++)
 167       for(m=0;(m<3);m++) {
 168         xcm[m]+=x[run_time][j][m]*w_rls[j];
 169       }
 170     for(m=0;(m<3);m++)
 171       xcm[m]/=tmas;
 172     for(j=0;(j<header.natoms);j++) {
 173       for(m=0;(m<3);m++)
 174         x[run_time][j][m]-=xcm[m];
 175     }
 176
 177     /*calculate energy of root_least_squares*/
 178     do_fit(header.natoms,w_rls,x[0],x[run_time]);
 179     fprintf(fp,"%8.4f %8.4f\n",
 180             t,calc_ener(header.natoms,w_rms,x[0],x[run_time]));
 181
 182     /*swap coordinates*/
 183     temp=x[0];
 184     for(i=0;(i<run_time);i++)
 185       x[i]=x[i+1];
 186     x[run_time]=temp;
 187   } while ((read_next_x(status,&t,natom,x[run_time],box)));
 188   fprintf(stderr,"\n");
 189   fclose(fp);
 190
 191   close_trj(status);
 192
 193   xvgr_file(ftp2fn(efXVG,NFILE,fnm),NULL);
 194
 195   thanx(stdout);
 196
 197   return 0;
 198 }