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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * GROtesk MACabre and Sinister
  31  */
  32
  33 #include <math.h>
  34 #include "smalloc.h"
  35 #include "typedefs.h"
  36 #include "copyrite.h"
  37 #include "statutil.h"
  38 #include "string2.h"
  39 #include "vec.h"
  40 #include "rdgroup.h"
  41 #include "pbc.h"
  42 #include "macros.h"
  43 #include "gstat.h"
  44 #include "xvgr.h"
  45 #include "tpxio.h"
  46
  47 real calc_ener(int natoms,real w_rms[],rvec x[],rvec xp[])
  48 {
  49   real e,tmas;
  50   int  j,m;
  51
  52   /*calculate energy */
  53   e=0;
  54   tmas=0;
  55   for(j=0;(j<natoms);j++) {
  56     tmas+=w_rms[j];
  57     for(m=0;(m<3);m++)
  58       e+=w_rms[j]*(x[j][m]-xp[j][m])*(x[j][m]-xp[j][m]);
  59   }
  60   /*return energy*/
  61   return(sqrt(e/tmas));
  62 }
  63
  64 int main (int argc,char *argv[])
  65 {
  66   static char *desc[] = {
  67     "g_run_rms computes the root mean square deviation (RMSD) of a structure",
  68     "from a trajectory x(t), with respect to a reference structure from the",
  69     "same trajectory, but at a specified time before the current time",
  70     "x(t-dt) by LSQ fitting the structures on top of each other[PAR]",
  71     "This tells you something about the mobility of structure as a function",
  72     "of time"
  73   };
  74   static int run_time = 5;
  75   t_pargs pa[] = {
  76     { "-t", FALSE, etINT, &run_time,
  77       "Time interval dt between reference and actual structure" }
  78   };
  79   int          step,nre,natom,i,j,m,teller=0;
  80   real         t,lambda,*w_rls,*w_rms,tmas;
  81   int          status;
  82   t_tpxheader  header;
  83   t_topology   top;
  84   matrix       box;
  85   rvec         **x,*xp,*v,xcm,*temp;
  86   FILE         *fp;
  87
  88 #define RLS 0
  89 #define RMS 1
  90   int          isize[2];
  91   atom_id      *index[2];
  92   char         *grpnames[2];
  93   t_filenm     fnm[] = {
  94     { efTRX, "-f", NULL, ffREAD },
  95     { efTPX, NULL, NULL, ffREAD },
  96     { efNDX, NULL, NULL, ffREAD },
  97     { efXVG, NULL, "runrms", ffWRITE }
  98   };
  99 #define NFILE asize(fnm)
 100
 101   CopyRight(stderr,argv[0]);
 102   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 103                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 104   snew(x,run_time+1);
 105
 106   read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&header,TRUE);
 107   snew(xp,header.natoms);
 108   for(i=0;(i<run_time+1);i++)
 109     snew(x[i],header.natoms);
 110   snew(w_rls,header.natoms);
 111   snew(w_rms,header.natoms);
 112   snew(v,header.natoms);
 113
 114   read_tpx(ftp2fn(efTPX,NFILE,fnm),
 115            &step,&t,&lambda,NULL,
 116            box,&natom,xp,NULL,NULL,&top);
 117
 118   /*set box type*/
 119   init_pbc(box,FALSE);
 120
 121   fprintf(stderr,"select group for root least squares fit\n");
 122   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize[RLS],&index[RLS],&grpnames[RLS]);
 123   for(i=0;(i<isize[RLS]);i++)
 124     w_rls[index[RLS][i]]=top.atoms.atom[index[RLS][i]].m;
 125
 126   fprintf(stderr,"select group for root mean square calculation\n");
 127   rd_index(ftp2fn(efNDX,NFILE,fnm),1,&isize[RMS],&index[RMS],&grpnames[RMS]);
 128   for(i=0;(i<isize[RMS]);i++)
 129     w_rms[index[RMS][i]]=top.atoms.atom[index[RMS][i]].m;
 130
 131   fp=xvgropen(ftp2fn(efXVG,NFILE,fnm),"Running RMSD","Time (ps)","RMSD (nm)");
 132
 133   rm_pbc(&(top.idef),top.atoms.nr,box,xp,xp);
 134
 135   /*calculate mass*/
 136   tmas=0;
 137   for(i=0;(i<header.natoms);i++)
 138     tmas+=w_rls[i];
 139
 140   /*set center of mass to zero for reference coordinates*/
 141
 142   /*do a first step*/
 143   natom=read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x[0],box);
 144   for(i=0; (i<run_time); i++) {
 145     read_next_x(status,&t,natom,x[i],box);
 146     rm_pbc(&(top.idef),top.atoms.nr,box,x[i],x[i]);
 147     for(m=0;(m<3);m++)
 148       xcm[m]=0;
 149     for(j=0;(j<header.natoms);j++)
 150       for(m=0;(m<3);m++) {
 151         xcm[m]+=x[i][j][m]*w_rls[j];
 152       }
 153     for(m=0;(m<3);m++)
 154       xcm[m]/=tmas;
 155     for(j=0;(j<header.natoms);j++) {
 156       for(m=0;(m<3);m++)
 157         x[i][j][m]-=xcm[m];
 158     }
 159   }
 160
 161   read_next_x(status,&t,natom,x[run_time],box);
 162   do {
 163     teller++;
 164     rm_pbc(&(top.idef),top.atoms.nr,box,x[run_time],x[run_time]);
 165     for(m=0;(m<3);m++)
 166       xcm[m]=0;
 167     for(j=0;(j<header.natoms);j++)
 168       for(m=0;(m<3);m++) {
 169         xcm[m]+=x[run_time][j][m]*w_rls[j];
 170       }
 171     for(m=0;(m<3);m++)
 172       xcm[m]/=tmas;
 173     for(j=0;(j<header.natoms);j++) {
 174       for(m=0;(m<3);m++)
 175         x[run_time][j][m]-=xcm[m];
 176     }
 177
 178     /*calculate energy of root_least_squares*/
 179     do_fit(header.natoms,w_rls,x[0],x[run_time]);
 180     fprintf(fp,"%8.4f %8.4f\n",
 181             t,calc_ener(header.natoms,w_rms,x[0],x[run_time]));
 182
 183     /*swap coordinates*/
 184     temp=x[0];
 185     for(i=0;(i<run_time);i++)
 186       x[i]=x[i+1];
 187     x[run_time]=temp;
 188   } while ((read_next_x(status,&t,natom,x[run_time],box)));
 189   fprintf(stderr,"\n");
 190   fclose(fp);
 191
 192   close_trj(status);
 193
 194   do_view(ftp2fn(efXVG,NFILE,fnm),NULL);
 195
 196   thanx(stderr);
 197
 198   return 0;
 199 }