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 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  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32
  33 #include "filenm.h"
  34 #include "smalloc.h"
  35 #include "macros.h"
  36 #include "math.h"
  37 #include "typedefs.h"
  38 #include "xvgr.h"
  39 #include "copyrite.h"
  40 #include "statutil.h"
  41 #include "tpxio.h"
  42 #include "string2.h"
  43 #include "vec.h"
  44 #include "rdgroup.h"
  45 #include "pbc.h"
  46 #include "fatal.h"
  47 #include "futil.h"
  48 #include "confio.h"
  49 #include "pdbio.h"
  50 #include "txtdump.h"
  51 #include "do_fit.h"
  52
  53 void rm_gropbc(t_atoms *atoms,rvec x[],matrix box)
  54 {
  55   real dist;
  56   int  n,d;
  57
  58   /* check periodic boundary */
  59   for(d=0; d<DIM; d++)
  60     for(n=1; n<atoms->nr; n++) {
  61       dist = x[n][d]-x[n-1][d];
  62       if ( fabs(dist) > 0.9 * box[d][d]  ) {
  63         if ( dist >  0 )
  64           x[n][d]-=box[d][d];
  65         else
  66           x[n][d]+=box[d][d];
  67       }
  68     }
  69 }
  70
  71 int main (int argc,char *argv[])
  72 {
  73   static char *desc[] = {
  74     "g_confrms computes the root mean square deviation (RMSD) of two",
  75     "structures after LSQ fitting the second structure on the first one.",
  76     "The two structures do NOT need to have the same number of atoms,",
  77     "only the two index groups used for the fit need to be identical.",
  78     "[PAR]",
  79     "The superimposed structures are written to file. In a [TT].pdb[tt] file",
  80     "the two structures will be written as separate models",
  81     "(use [TT]rasmol -nmrpdb[tt]).",
  82   };
  83   static bool bOne=FALSE,bRmpbc=FALSE,bFit=TRUE;
  84
  85   t_pargs pa[] = {
  86     { "-one", FALSE, etBOOL, {&bOne},   "Only write the fitted structure to file" },
  87     { "-pbc", FALSE, etBOOL, {&bRmpbc}, "Try to make molecules whole again" },
  88     { "-fit", FALSE, etBOOL, {&bFit},
  89       "Do least squares superposition of the target structure to the reference" },
  90   };
  91   t_filenm fnm[] = {
  92     { efTPS, "-f1",  "conf1.gro", ffREAD  },
  93     { efSTX, "-f2",  "conf2",     ffREAD  },
  94     { efSTO, "-o",   "fit.pdb",   ffWRITE },
  95     { efNDX, "-n1" , "fit1.ndx",  ffOPTRD },
  96     { efNDX, "-n2" , "fit2.ndx",  ffOPTRD }
  97   };
  98 #define NFILE asize(fnm)
  99
 100   /* the two structure files */
 101   FILE    *fp;
 102   char    title1[STRLEN],title2[STRLEN],*name1,*name2;
 103   t_topology top;
 104   t_atoms atoms1,atoms2;
 105   int     natoms1,natoms2,warn=0;
 106   matrix  box;
 107   atom_id at;
 108   real    *w_rls,mass,totmass;
 109   rvec    *x1,*v1,*x2,*v2,*fit_x;
 110   matrix  box1,box2;
 111
 112   /* counters */
 113   int     i,j,m;
 114
 115   /* center of mass calculation */
 116   real    tmas1,tmas2;
 117   rvec    xcm1,xcm2;
 118
 119   /* variables for fit */
 120   char    *groupnames1,*groupnames2;
 121   int     isize1,isize2;
 122   atom_id *index1,*index2;
 123   real    rms;
 124
 125   CopyRight(stderr,argv[0]);
 126   parse_common_args(&argc,argv,PCA_BE_NICE,
 127                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 128
 129   /* reading reference structure from first structure file */
 130   fprintf(stderr,"\nReading first structure file\n");
 131   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title1,&top,&x1,&v1,box,TRUE);
 132   atoms1 = top.atoms;
 133   fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
 134           title1,atoms1.nr,atoms1.nres);
 135   srenew(atoms1.resname,atoms1.nres);
 136
 137   if ( bRmpbc )
 138     rm_gropbc(&atoms1,x1,box1);
 139
 140   fprintf(stderr,"Select group from first structure\n");
 141   get_index(&atoms1,opt2fn_null("-n1",NFILE,fnm),
 142             1,&isize1,&index1,&groupnames1);
 143
 144   if (bFit && (isize1 < 3))
 145     fatal_error(0,"Need >= 3 points to fit!\n");
 146
 147   /* reading second structure file */
 148   get_stx_coordnum(opt2fn("-f2",NFILE,fnm),&(atoms2.nr));
 149   snew(x2,atoms2.nr);
 150   snew(v2,atoms2.nr);
 151   snew(atoms2.resname,atoms2.nr);
 152   snew(atoms2.atom,atoms2.nr);
 153   snew(atoms2.atomname,atoms2.nr);
 154   fprintf(stderr,"\nReading second structure file\n");
 155   read_stx_conf(opt2fn("-f2",NFILE,fnm),title2,&atoms2,x2,v2,box2);
 156   fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
 157           title2,atoms2.nr,atoms2.nres);
 158   srenew(atoms2.resname,atoms2.nres);
 159
 160   if ( bRmpbc )
 161     rm_gropbc(&atoms2,x2,box2);
 162
 163   fprintf(stderr,"Select group from second structure\n");
 164   get_index(&atoms2,opt2fn_null("-n2",NFILE,fnm),
 165             1,&isize2,&index2,&groupnames2);
 166
 167   /* check isizes, must be equal */
 168   if ( isize2 != isize1 )
 169     fatal_error(0,"You selected groups with differen number of atoms.\n");
 170
 171   for(i=0; i<isize1; i++) {
 172     name1=*atoms1.atomname[index1[i]];
 173     name2=*atoms2.atomname[index2[i]];
 174     if (strcmp(name1,name2)) {
 175       if (warn < 20)
 176         fprintf(stderr,
 177                 "Warning: atomnames at index %d don't match: %u %s, %u %s\n",
 178                 i+1,index1[i]+1,name1,index2[i]+1,name2);
 179       warn++;
 180     }
 181   }
 182   if (warn)
 183     fprintf(stderr,"%d atomname%s did not match\n",warn,(warn==1) ? "":"s");
 184
 185   if (bFit) {
 186     /* calculate center of geometry of reference structure */
 187     for(m=0; m<DIM; m++) {
 188       xcm1[m]=0;
 189       for(i=0; i<isize1; i++)
 190         xcm1[m]+=x1[index1[i]][m];
 191       xcm1[m]/=isize1;
 192       for(i=0; i<atoms1.nr; i++)
 193         x1[i][m]-=xcm1[m];
 194     }
 195
 196     /* calculate center of geometry of structure to be fitted */
 197     for(m=0; m<DIM; m++) {
 198       xcm2[m]=0;
 199       for(i=0; i<isize2; i++)
 200         xcm2[m]+=x2[index2[i]][m];
 201       xcm2[m]/=isize2;
 202       for(i=0; i<atoms2.nr; i++)
 203         x2[i][m]-=xcm2[m];
 204     }
 205
 206     snew(w_rls,atoms2.nr);
 207     snew(fit_x,atoms2.nr);
 208     for(at=0; at<isize1; at++) {
 209       w_rls[index2[at]] = 1.0;
 210       copy_rvec(x1[index1[at]],fit_x[index2[at]]);
 211     }
 212
 213     /* do the least squares fit to the reference structure */
 214     do_fit(atoms2.nr,w_rls,fit_x,x2);
 215
 216     sfree(fit_x);
 217     sfree(w_rls);
 218   }
 219   else {
 220     clear_rvec(xcm1);
 221     clear_rvec(xcm2);
 222   }
 223
 224   /* calculate the rms deviation */
 225   rms     = 0;
 226   for(at=0; at<isize1; at++) {
 227     for(m=0; m<DIM; m++)
 228       rms += sqr(x1[index1[at]][m] - x2[index2[at]][m]);
 229   }
 230   rms = sqrt(rms/isize1);
 231
 232   fprintf(stderr,"Root mean square deviation after lsq fit = %g\n",rms);
 233
 234   if (bFit) {
 235     /* reset coordinates of reference and fitted structure */
 236     for(i=0; i<atoms1.nr; i++)
 237       for(m=0; m<DIM; m++)
 238         x1[i][m]+=xcm1[m];
 239     for(i=0; i<atoms2.nr; i++)
 240       for(m=0; m<DIM; m++)
 241         x2[i][m]+=xcm1[m];
 242   }
 243   /* write gromos file of fitted structure(s) */
 244   fp=ffopen(opt2fn("-o",NFILE,fnm),"w");
 245   if (fn2ftp(opt2fn("-o",NFILE,fnm))==efGRO) {
 246     if (!bOne)
 247       write_hconf_p(fp,title1,&atoms1,3,x1,v1,box1);
 248     write_hconf_p(fp,title2,&atoms2,3,x2,v2,box2);
 249   } else {
 250     if (!bOne)
 251       write_pdbfile(fp,title1,&atoms1,x1,box1,0,1);
 252     write_pdbfile(fp,title2,&atoms2,x2,box2,0,bOne ? -1 : 2);
 253   }
 254   fclose(fp);
 255
 256   thanx(stderr);
 257
 258   return 0;
 259 }