src/tools/g_confrms.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_confrms_c = "$Id$";
  30
  31 /*
  32  * fit coordinates of file 1 to coordinates of file 2
  33  * fit is done on the atoms as indicated by the NDX- file
  34  */
  35
  36 #include "filenm.h"
  37 #include "smalloc.h"
  38 #include "macros.h"
  39 #include "math.h"
  40 #include "typedefs.h"
  41 #include "xvgr.h"
  42 #include "copyrite.h"
  43 #include "statutil.h"
  44 #include "tpxio.h"
  45 #include "string2.h"
  46 #include "vec.h"
  47 #include "rdgroup.h"
  48 #include "pbc.h"
  49 #include "fatal.h"
  50 #include "futil.h"
  51 #include "confio.h"
  52 #include "pdbio.h"
  53 #include "txtdump.h"
  54 #include "do_fit.h"
  55
  56 #define EPS  1.0e-09
  57
  58 void rm_gropbc(t_atoms *atoms,rvec x[],matrix box)
  59 {
  60   real dist;
  61   int  n,d;
  62
  63   /* check periodic boundary */
  64   for(d=0;(d<DIM);d++) {
  65     for(n=1;(n<atoms->nr);n++) {
  66       dist = x[n][d]-x[n-1][d];
  67       if ( fabs(dist) > 0.9 * box[d][d]  ) {
  68         if ( dist >  0 )
  69           x[n][d]-=box[d][d];
  70         else
  71           x[n][d]+=box[d][d];
  72       }
  73     }
  74   }
  75 }
  76
  77 int main (int argc,char *argv[])
  78 {
  79   static char *desc[] = {
  80     "g_confrms computes the root mean square deviation (RMSD) of two",
  81     "structures after LSQ fitting the second structure on the first one.",
  82     "The two structures do NOT need to have the same number of atoms,",
  83     "only the two index groups used for the fit need to be identical.",
  84     "[PAR]",
  85     "The superimposed structures are written to file. In a [TT].pdb[tt] file",
  86     "the two structures will have chain identifiers 'A' and 'B' respectively.",
  87     "When the option [TT]-one[tt] is set, only the fitted structure is",
  88     "written to file and the chain identifiers are not changed."
  89   };
  90   static bool bSecond=FALSE,bRmpbc=FALSE;
  91
  92   t_pargs pa[] = {
  93     { "-one", FALSE, etBOOL, {&bSecond}, "Only write the fitted structure to file" },
  94     { "-pbc", FALSE, etBOOL, {&bRmpbc}, "Try to make molecules whole again" }
  95   };
  96   t_filenm fnm[] = {
  97     { efTPS, "-f1",  "conf1.gro", ffREAD  },
  98     { efSTX, "-f2",  "conf2", ffREAD  },
  99     { efSTO, "-o",    "fit.pdb",  ffWRITE },
 100     { efNDX, "-n1" ,  "fit1.ndx",   ffOPTRD  },
 101     { efNDX, "-n2" ,  "fit2.ndx",   ffOPTRD  }
 102   };
 103
 104 #define NFILE asize(fnm)
 105
 106   /* the two gromos files */
 107   FILE    *fp;
 108   char    title_1[STRLEN],title_2[STRLEN],*name1,*name2;
 109   t_topology top;
 110   t_atoms atoms_1,atoms_2;
 111   int     natoms_1,natoms_2,warn=0;
 112   matrix  box;
 113   atom_id at;
 114   real    *w_rls,mass,totmass;
 115   rvec    *x_1,*v_1,*x_2,*v_2,*xf_1,*xf_2,*fit_x;
 116   matrix  box_1,box_2;
 117
 118   /* counters */
 119   int   i,j,m;
 120
 121   /* center of mass calculation */
 122   real tmas_1,tmas_2;
 123   rvec xcm_1,xcm_2;
 124
 125   /* variables for fit */
 126   char *groupnames_1,*groupnames_2;
 127   int isize_1,isize_2;
 128   atom_id *index_1,*index_2;
 129   real rms;
 130
 131   CopyRight(stderr,argv[0]);
 132   parse_common_args(&argc,argv,0,TRUE,
 133                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 134
 135   /* reading first structure file, nb this is going
 136    * to be the reference structure*/
 137   fprintf(stderr,"\nReading first structure file\n");
 138   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title_1,&top,&x_1,&v_1,box,TRUE);
 139   atoms_1 = top.atoms;
 140   fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
 141           title_1,atoms_1.nr,atoms_1.nres);
 142   srenew(atoms_1.resname,atoms_1.nres);
 143
 144   if ( bRmpbc )
 145     rm_gropbc(&atoms_1,x_1,box_1);
 146
 147
 148   /* check if we have an NDX file */
 149   /*
 150   if ( opt2bSet("-n1",NFILE,fnm) ) {
 151     fprintf(stderr,"\nSelect group for root least square fit\n");
 152     rd_index(opt2fn("-n1",NFILE,fnm),1,&isize_1,&index_1,&groupnames_1);
 153   } else {
 154     isize_1 = atoms_1.nr;
 155     snew(index_1,isize_1);
 156     for(i=0;(i<isize_1);i++)
 157       index_1[i]=i;
 158     groupnames_1 = title_1;
 159   }
 160   */
 161   get_index(&atoms_1,opt2fn_null("-n1",NFILE,fnm),
 162             1,&isize_1,&index_1,&groupnames_1);
 163
 164   if (isize_1 < 3)
 165     fatal_error(0,"Need >= 3 points to fit!\n");
 166
 167
 168   /* reading second gromos file */
 169   get_stx_coordnum(opt2fn("-f2",NFILE,fnm),&(atoms_2.nr));
 170   snew(x_2,atoms_2.nr);
 171   snew(v_2,atoms_2.nr);
 172   snew(atoms_2.resname,atoms_2.nr);
 173   snew(atoms_2.atom,atoms_2.nr);
 174   snew(atoms_2.atomname,atoms_2.nr);
 175   fprintf(stderr,"\nReading second structure file\n");
 176   read_stx_conf(opt2fn("-f2",NFILE,fnm),title_2,&atoms_2,x_2,v_2,box_2);
 177   fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
 178           title_2,atoms_2.nr,atoms_2.nres);
 179   srenew(atoms_2.resname,atoms_2.nres);
 180
 181   if ( bRmpbc )
 182     rm_gropbc(&atoms_2,x_2,box_2);
 183
 184   get_index(&atoms_2,opt2fn_null("-n2",NFILE,fnm),
 185             1,&isize_2,&index_2,&groupnames_2);
 186
 187   /* check isizes, must be equal */
 188   if ( isize_2 != isize_1 )
 189     fatal_error(0,"isize_2 != isize_2\n");
 190
 191   if (isize_2 < 3)
 192     fatal_error(0,"Need >= 3 points to fit!\n");
 193
 194   for(i=0; (i<isize_1); i++) {
 195     name1=*atoms_1.atomname[index_1[i]];
 196     name2=*atoms_2.atomname[index_2[i]];
 197     if (strcmp(name1,name2)) {
 198       if (warn<20)
 199         fprintf(stderr,
 200                 "Warning: atomnames at index %d don't match: %u %s, %u %s\n",
 201                 i+1,index_1[i]+1,name1,index_2[i]+1,name2);
 202       warn++;
 203     }
 204   }
 205   if (warn)
 206     fprintf(stderr,"%d atomname%s did not match\n",warn,(warn==1) ? "":"s");
 207
 208   /* calculate center of mass of reference structure */
 209   for(m=0;(m<3);m++)
 210     xcm_1[m]=0;
 211   for(i=0;(i<isize_1);i++)
 212     for(m=0;(m<DIM);m++)
 213       xcm_1[m]+=x_1[index_1[i]][m];
 214   for(m=0;(m<3);m++)
 215     xcm_1[m]/=isize_1;
 216   for(i=0;(i<atoms_1.nr);i++)
 217     for(m=0;(m<DIM);m++)
 218       x_1[i][m]-=xcm_1[m];
 219
 220   /* calculate center of mass of structure to be fitted */
 221   for(m=0;(m<3);m++)
 222     xcm_2[m]=0;
 223   for(i=0;(i<isize_2);i++)
 224     for(m=0;(m<DIM);m++)
 225       xcm_2[m]+=x_2[index_2[i]][m];
 226   for(m=0;(m<3);m++)
 227     xcm_2[m]/=isize_2;
 228   for(i=0;(i<atoms_2.nr);i++)
 229     for(m=0;(m<DIM);m++)
 230       x_2[i][m]-=xcm_2[m];
 231
 232   snew(w_rls,atoms_2.nr);
 233   snew(fit_x,atoms_2.nr);
 234   for(at=0; at<isize_1; at++) {
 235     w_rls[index_2[at]] = 1.0;
 236     copy_rvec(x_1[index_1[at]],fit_x[index_2[at]]);
 237   }
 238
 239   /*do the least squares fit to the reference structure*/
 240   /*
 241   rms = my_fit(isize_1,index_1,index_2,x_1,atoms_2.nr,x_2);
 242   */
 243   do_fit(atoms_2.nr,w_rls,fit_x,x_2);
 244
 245   rms     = 0;
 246   totmass = 0;
 247   for(at=0; at<isize_1; at++) {
 248     mass     = w_rls[index_2[at]];
 249     for(m=0; m<3; m++)
 250       rms     += sqr(x_1[index_1[at]][m] - x_2[index_2[at]][m])*mass;
 251     totmass += mass;
 252   }
 253   rms = sqrt(rms/totmass);
 254
 255   fprintf(stderr,"Root mean square deviation after lsq fit = %g\n",rms);
 256
 257   /* reset coordinates of reference and fitted structure */
 258   for(i=0;(i<atoms_1.nr);i++) {
 259     for(m=0;(m<3);m++)
 260       x_1[i][m]+=xcm_1[m];
 261   }
 262   for(i=0;(i<atoms_2.nr);i++) {
 263     for(m=0;(m<3);m++)
 264       x_2[i][m]+=xcm_1[m];
 265   }
 266
 267   /* calculate the rms deviation */
 268
 269
 270
 271   /* write gromos file of fitted structure(s) */
 272   fp=ffopen(opt2fn("-o",NFILE,fnm),"w");
 273   if (fn2ftp(opt2fn("-o",NFILE,fnm))==efGRO) {
 274     if (!bSecond)
 275       write_hconf(fp,title_1,&atoms_1,x_1,v_1,box_1);
 276     write_hconf(fp,title_2,&atoms_2,x_2,v_2,box_2);
 277   } else {
 278     if (bSecond)
 279       write_pdbfile(fp,title_2,&atoms_2,x_2,box_2,0,TRUE);
 280     else {
 281       write_pdbfile(fp,title_1,&atoms_1,x_1,box_1,'A',FALSE);
 282       write_pdbfile(fp,title_2,&atoms_2,x_2,box_2,'B',TRUE);
 283     }
 284   }
 285   fclose(fp);
 286
 287   thanx(stdout);
 288
 289   return 0;
 290 }
 291
 292
 293
 294
 295
 296
 297
 298
 299