src/tools/gmx_confrms.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "filenm.h"
  43 #include "smalloc.h"
  44 #include "macros.h"
  45 #include <math.h>
  46 #include "typedefs.h"
  47 #include "xvgr.h"
  48 #include "copyrite.h"
  49 #include "statutil.h"
  50 #include "tpxio.h"
  51 #include "string2.h"
  52 #include "vec.h"
  53 #include "index.h"
  54 #include "pbc.h"
  55 #include "gmx_fatal.h"
  56 #include "futil.h"
  57 #include "confio.h"
  58 #include "pdbio.h"
  59 #include "txtdump.h"
  60 #include "do_fit.h"
  61 #include "viewit.h"
  62 #include "rmpbc.h"
  63 #include "gmx_ana.h"
  64
  65
  66 void calc_rm_cm(int isize, atom_id index[], t_atoms *atoms, rvec x[], rvec xcm)
  67 {
  68   int i,d;
  69   real tm,m;
  70
  71   /* calculate and remove center of mass of reference structure */
  72   tm = 0;
  73   clear_rvec(xcm);
  74   for(i=0; i<isize; i++) {
  75     m = atoms->atom[index[i]].m;
  76     for(d=0; d<DIM; d++)
  77       xcm[d] += m*x[index[i]][d];
  78     tm += m;
  79   }
  80   svmul(1/tm,xcm,xcm);
  81   for(i=0; i<atoms->nr; i++)
  82     rvec_dec(x[i],xcm);
  83 }
  84
  85 int build_res_index(int isize, atom_id index[], t_atom atom[], int rindex[])
  86 {
  87   int i, r;
  88
  89   r=0;
  90   rindex[r] = atom[index[0]].resind;
  91   r++;
  92   for(i=1; i<isize; i++)
  93     if (atom[index[i]].resind != rindex[r-1]) {
  94       rindex[r] = atom[index[i]].resind;
  95       r++;
  96     }
  97
  98   return r;
  99 }
 100
 101 int find_res_end(int i, int isize, atom_id index[], t_atoms *atoms)
 102 {
 103   int rnr;
 104
 105   rnr = atoms->atom[index[i]].resind;
 106   while(i<isize && atoms->atom[index[i]].resind==rnr)
 107     i++;
 108   return i;
 109 }
 110
 111 int debug_strcmp(char s1[], char s2[])
 112 {
 113   if (debug) fprintf(debug," %s-%s",s1,s2);
 114   return strcmp(s1, s2);
 115 }
 116
 117 int find_next_match_atoms_in_res(int *i1,int isize1,atom_id index1[],
 118                                  int m1,char **atnms1[],
 119                                  int *i2,int isize2,atom_id index2[],
 120                                  int m2,char **atnms2[])
 121 {
 122   int dx, dy, dmax, cmp;
 123   gmx_bool bFW=FALSE;
 124
 125   dx=dy=0;
 126   cmp=NOTSET;
 127   dmax = max(m1-*i1, m2-*i2);
 128   for(dx=0; dx<dmax && cmp!=0; dx++) {
 129     for(dy=dx; dy<dmax && cmp!=0; dy++) {
 130       if (dx || dy) {
 131         if (debug) fprintf(debug, ".");
 132         cmp=NOTSET;
 133         if ( *i1+dx<m1 && *i2+dy<m2 ) {
 134           bFW=TRUE;
 135           cmp = debug_strcmp(*atnms1[index1[*i1+dx]],*atnms2[index2[*i2+dy]]);
 136           if (debug) fprintf(debug,"(%d %d)", *i1+dx, *i2+dy);
 137         }
 138         if ( cmp!=0 && *i1+dy<m1 && *i2+dx<m2 ) {
 139           bFW=FALSE;
 140           cmp = debug_strcmp(*atnms1[index1[*i1+dy]],*atnms2[index2[*i2+dx]]);
 141           if (debug) fprintf(debug,"(%d %d)", *i1+dy, *i2+dx);
 142         }
 143       }
 144     }
 145   }
 146   /* apparently, dx and dy are incremented one more time
 147      as the loops terminate: we correct this here */
 148   dx--;
 149   dy--;
 150   if (cmp==0) {
 151     if (debug) fprintf(debug, "{%d %d}", *i1+bFW?dx:dy, *i2+bFW?dy:dx );
 152     if (bFW) {
 153       *i1 += dx;
 154       *i2 += dy;
 155     } else {
 156       *i1 += dy;
 157       *i2 += dx;
 158     }
 159   }
 160
 161   return cmp;
 162 }
 163
 164 static int find_next_match_res(int *rnr1, int isize1,
 165                                int index1[], t_resinfo *resinfo1,
 166                                int *rnr2, int isize2,
 167                                int index2[], t_resinfo *resinfo2)
 168 {
 169   int dx, dy, dmax, cmp, rr1, rr2;
 170   gmx_bool bFW=FALSE,bFF=FALSE;
 171
 172   dx=dy=0;
 173   rr1 = 0;
 174   while(rr1<isize1 && *rnr1 != index1[rr1])
 175     rr1++;
 176   rr2 = 0;
 177   while(rr2<isize2 && *rnr2 != index2[rr2])
 178     rr2++;
 179
 180   cmp=NOTSET;
 181   dmax = max(isize1-rr1, isize2-rr2);
 182   if (debug) fprintf(debug," R:%d-%d:%d-%d:%d ",
 183                      rr1, isize1, rr2, isize2, dmax);
 184   for(dx=0; dx<dmax && cmp!=0; dx++)
 185     for(dy=0; dy<=dx && cmp!=0; dy++)
 186       if ( dx!=dy ) {
 187         cmp=NOTSET;
 188         if ( rr1+dx<isize1 && rr2+dy<isize2 ) {
 189           bFW=TRUE;
 190           cmp=debug_strcmp(*resinfo1[index1[rr1+dx]].name,
 191                            *resinfo2[index2[rr2+dy]].name);
 192           if (debug) fprintf(debug,"(%d %d)", rr1+dx, rr2+dy);
 193         }
 194         if ( cmp!=0 && rr1+dy<isize1 && rr2+dx<isize2 ) {
 195           bFW=FALSE;
 196           cmp=debug_strcmp(*resinfo1[index1[rr1+dy]].name,
 197                            *resinfo2[index2[rr2+dx]].name);
 198           if (debug) fprintf(debug,"(%d %d)", rr1+dy, rr2+dx);
 199         }
 200         if ( dx!=0 && cmp!=0 && rr1+dx<isize1 && rr2+dx<isize2 ) {
 201           bFF=TRUE;
 202           cmp=debug_strcmp(*resinfo1[index1[rr1+dx]].name,
 203                            *resinfo2[index2[rr2+dx]].name);
 204           if (debug) fprintf(debug,"(%d %d)", rr1+dx, rr2+dx);
 205         } else
 206           bFF=FALSE;
 207       }
 208   /* apparently, dx and dy are incremented one more time
 209      as the loops terminate: we correct this here */
 210   dx--;
 211   dy--;
 212   /* if we skipped equal on both sides, only skip one residue
 213      to allow for single mutations of residues... */
 214   if (bFF) {
 215     if (debug) fprintf(debug, "%d.%d.%dX%sX%s",dx,rr1,rr2,
 216                        *resinfo1[index1[rr1+1]].name,
 217                        *resinfo2[index2[rr2+1]].name);
 218     dx=1;
 219   }
 220   if (cmp==0) {
 221     if (debug) fprintf(debug, "!");
 222     if (bFF) {
 223       rr1 += dx;
 224       rr2 += dx;
 225     } else
 226       if (bFW) {
 227         rr1 += dx;
 228         rr2 += dy;
 229       } else {
 230         rr1 += dy;
 231         rr2 += dx;
 232       }
 233     *rnr1 = index1[rr1];
 234     *rnr2 = index2[rr2];
 235   }
 236
 237   return cmp;
 238 }
 239
 240 int find_first_atom_in_res(int rnr, int isize, atom_id index[], t_atom atom[])
 241 {
 242   int i;
 243
 244   i=0;
 245   while(i<isize && atom[index[i]].resind!=rnr)
 246     i++;
 247
 248   if (atom[index[i]].resind==rnr)
 249     return i;
 250   else
 251     return NOTSET;
 252 }
 253
 254 void find_matching_names(int *isize1, atom_id index1[], t_atoms *atoms1,
 255                          int *isize2, atom_id index2[], t_atoms *atoms2)
 256 {
 257   int i, i1, i2, ii1, ii2, m1, m2;
 258   int atcmp, rescmp;
 259   int r,rnr1, rnr2, prnr1, prnr2;
 260   int rsize1, rsize2;
 261   int *rindex1, *rindex2;
 262   char *resnm1, *resnm2, *atnm1, *atnm2;
 263   char ***atnms1,***atnms2;
 264   t_resinfo *resinfo1,*resinfo2;
 265
 266   /* set some handy shortcuts */
 267   resinfo1 = atoms1->resinfo;
 268   atnms1   = atoms1->atomname;
 269   resinfo2 = atoms2->resinfo;
 270   atnms2   = atoms2->atomname;
 271
 272   /* build indexes of selected residues */
 273   snew(rindex1,atoms1->nres);
 274   rsize1=build_res_index(*isize1, index1, atoms1->atom, rindex1);
 275   snew(rindex2,atoms2->nres);
 276   rsize2=build_res_index(*isize2, index2, atoms2->atom, rindex2);
 277
 278   i1=i2=0;
 279   ii1=ii2=0;
 280   atcmp=rescmp=0;
 281   prnr1=prnr2=NOTSET;
 282   if (debug) fprintf(debug, "Find matching names: %d, %d\n", *isize1, *isize2);
 283   while ( atcmp==0 && i1<*isize1 && i2<*isize2 ) {
 284     /* prologue */
 285     rnr1 = atoms1->atom[index1[i1]].resind;
 286     rnr2 = atoms2->atom[index2[i2]].resind;
 287     if (rnr1!=prnr1 || rnr2!=prnr2) {
 288       if (debug) fprintf(debug, "R: %s%d %s%d\n",
 289                          *resinfo1[rnr1].name,rnr1,*resinfo2[rnr2].name,rnr2);
 290       rescmp=strcmp(*resinfo1[rnr1].name, *resinfo2[rnr2].name);
 291     }
 292     if (debug) fprintf(debug, "comparing %d %d", i1, i2);
 293     atcmp=debug_strcmp(*atnms1[index1[i1]], *atnms2[index2[i2]]);
 294
 295     /* the works */
 296     if (atcmp!=0) { /* no match -> find match within residues */
 297       m1 = find_res_end(i1, *isize1, index1, atoms1);
 298       m2 = find_res_end(i2, *isize2, index2, atoms2);
 299       if (debug) fprintf(debug, " [%d<%d %d<%d]", i1,m1, i2,m2);
 300       atcmp = find_next_match_atoms_in_res(&i1,*isize1,index1,m1,atnms1,
 301                                            &i2,*isize2,index2,m2,atnms2);
 302       if (debug) fprintf(debug, " -> %d %d %s-%s", i1, i2,
 303                          *atnms1[index1[i1]],*atnms2[index2[i2]]);
 304
 305     }
 306     if (atcmp!=0) { /* still no match -> next residue(s) */
 307       prnr1=rnr1;
 308       prnr2=rnr2;
 309       rescmp = find_next_match_res(&rnr1,rsize1,rindex1,resinfo1,
 310                                    &rnr2,rsize2,rindex2,resinfo2);
 311       if (rnr1!=prnr1)
 312         i1 = find_first_atom_in_res(rnr1, *isize1, index1, atoms1->atom);
 313       if (rnr2!=prnr2)
 314         i2 = find_first_atom_in_res(rnr2, *isize2, index2, atoms2->atom);
 315       if (debug) fprintf(debug, " -> %s%d-%s%d %s%d-%s%d",
 316                          *resinfo1[rnr1].name,rnr1,
 317                          *atnms1[index1[i1]],index1[i1],
 318                          *resinfo2[rnr2].name,rnr2,
 319                          *atnms2[index2[i2]],index2[i2]);
 320       m1 = find_res_end(i1, *isize1, index1, atoms1);
 321       m2 = find_res_end(i2, *isize2, index2, atoms2);
 322       if (debug) fprintf(debug, " [%d<%d %d<%d]", i1,m1, i2,m2);
 323       atcmp = find_next_match_atoms_in_res(&i1,*isize1,index1,m1,atnms1,
 324                                            &i2,*isize2,index2,m2,atnms2);
 325       if (debug) fprintf(debug, " -> %d %d %s-%s", i1, i2,
 326                          *atnms1[index1[i1]],*atnms2[index2[i2]]);
 327     }
 328     if (debug)
 329       fprintf(debug, "(%d %d): %d %d\n", ii1, ii2, atcmp, rescmp);
 330     if (atcmp==0) { /* if match -> store indices */
 331       index1[ii1++]=index1[i1];
 332       index2[ii2++]=index2[i2];
 333     }
 334     i1++;
 335     i2++;
 336
 337     /* epilogue */
 338     prnr1=rnr1;
 339     prnr2=rnr2;
 340   }
 341
 342   if (ii1 != ii2)
 343     gmx_fatal(FARGS, "DEATH HORROR: non-equal number of matching atoms!\n");
 344   if ( ii1==i1 && ii2==i2 )
 345     printf("All atoms in index groups 1 and 2 match\n");
 346   else {
 347     if ( i1==i2 && ii1==ii2 )
 348       printf("Both index groups modified from %d to %d atoms\n", i1, ii1);
 349     else {
 350       if (ii1 != i1)
 351         printf("Index group 1 modified from %d to %d atoms\n",i1, ii1);
 352       if (ii2 != i2)
 353         printf("Index group 2 modified from %d to %d atoms\n",i2, ii2);
 354     }
 355     *isize1 = ii1;
 356     *isize2 = ii2;
 357   }
 358 }
 359 /* 1 */
 360
 361 int gmx_confrms(int argc,char *argv[])
 362 {
 363   const char *desc[] = {
 364     "[TT]g_confrms[tt] computes the root mean square deviation (RMSD) of two",
 365     "structures after least-squares fitting the second structure on the first one.",
 366     "The two structures do NOT need to have the same number of atoms,",
 367     "only the two index groups used for the fit need to be identical.",
 368     "With [TT]-name[tt] only matching atom names from the selected groups",
 369     "will be used for the fit and RMSD calculation. This can be useful ",
 370     "when comparing mutants of a protein.",
 371     "[PAR]",
 372     "The superimposed structures are written to file. In a [TT].pdb[tt] file",
 373     "the two structures will be written as separate models",
 374     "(use [TT]rasmol -nmrpdb[tt]). Also in a [TT].pdb[tt] file, B-factors",
 375     "calculated from the atomic MSD values can be written with [TT]-bfac[tt].",
 376   };
 377   static gmx_bool bOne=FALSE,bRmpbc=FALSE,bMW=TRUE,bName=FALSE,
 378     bBfac=FALSE,bFit=TRUE,bLabel=FALSE;
 379
 380   t_pargs pa[] = {
 381     { "-one", FALSE, etBOOL, {&bOne},   "Only write the fitted structure to file" },
 382     { "-mw",  FALSE, etBOOL, {&bMW},    "Mass-weighted fitting and RMSD" },
 383     { "-pbc", FALSE, etBOOL, {&bRmpbc}, "Try to make molecules whole again" },
 384     { "-fit", FALSE, etBOOL, {&bFit},
 385       "Do least squares superposition of the target structure to the reference" },
 386     { "-name", FALSE,etBOOL,{&bName},
 387       "Only compare matching atom names" },
 388     { "-label",FALSE,etBOOL,{&bLabel},
 389       "Added chain labels A for first and B for second structure"},
 390     { "-bfac", FALSE,etBOOL,{&bBfac},
 391       "Output B-factors from atomic MSD values" }
 392   };
 393   t_filenm fnm[] = {
 394     { efTPS, "-f1",  "conf1.gro", ffREAD  },
 395     { efSTX, "-f2",  "conf2",     ffREAD  },
 396     { efSTO, "-o",   "fit.pdb",   ffWRITE },
 397     { efNDX, "-n1" , "fit1.ndx",  ffOPTRD },
 398     { efNDX, "-n2", "fit2.ndx",  ffOPTRD },
 399     { efNDX, "-no", "match.ndx", ffOPTWR }
 400   };
 401 #define NFILE asize(fnm)
 402
 403   /* the two structure files */
 404   const char *conf1file, *conf2file, *matchndxfile, *outfile;
 405   FILE    *fp;
 406   char    title1[STRLEN],title2[STRLEN],*name1,*name2;
 407   t_topology *top1,*top2;
 408   int     ePBC1,ePBC2;
 409   t_atoms *atoms1,*atoms2;
 410   int     warn=0;
 411   atom_id at;
 412   real    *w_rls,mass,totmass;
 413   rvec    *x1,*v1,*x2,*v2,*fit_x;
 414   matrix  box1,box2;
 415
 416   output_env_t oenv;
 417
 418   /* counters */
 419   int     i,j,m;
 420
 421   /* center of mass calculation */
 422   real    tmas1,tmas2;
 423   rvec    xcm1,xcm2;
 424
 425   /* variables for fit */
 426   char    *groupnames1,*groupnames2;
 427   int     isize1,isize2;
 428   atom_id *index1,*index2;
 429   real    rms,msd,minmsd,maxmsd;
 430   real    *msds;
 431
 432
 433   CopyRight(stderr,argv[0]);
 434   parse_common_args(&argc,argv,PCA_BE_NICE | PCA_CAN_VIEW,
 435                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 436   matchndxfile = opt2fn_null("-no",NFILE,fnm);
 437   conf1file = ftp2fn(efTPS,NFILE,fnm);
 438   conf2file = ftp2fn(efSTX,NFILE,fnm);
 439
 440   /* reading reference structure from first structure file */
 441   fprintf(stderr,"\nReading first structure file\n");
 442   snew(top1,1);
 443   read_tps_conf(conf1file,title1,top1,&ePBC1,&x1,&v1,box1,TRUE);
 444   atoms1 = &(top1->atoms);
 445   fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
 446           title1,atoms1->nr,atoms1->nres);
 447
 448   if ( bRmpbc )
 449     rm_gropbc(atoms1,x1,box1);
 450
 451   fprintf(stderr,"Select group from first structure\n");
 452   get_index(atoms1,opt2fn_null("-n1",NFILE,fnm),
 453             1,&isize1,&index1,&groupnames1);
 454   printf("\n");
 455
 456   if (bFit && (isize1 < 3))
 457     gmx_fatal(FARGS,"Need >= 3 points to fit!\n");
 458
 459   /* reading second structure file */
 460   fprintf(stderr,"\nReading second structure file\n");
 461   snew(top2,1);
 462   read_tps_conf(conf2file,title2,top2,&ePBC2,&x2,&v2,box2,TRUE);
 463   atoms2 = &(top2->atoms);
 464   fprintf(stderr,"%s\nContaining %d atoms in %d residues\n",
 465           title2,atoms2->nr,atoms2->nres);
 466
 467   if ( bRmpbc )
 468     rm_gropbc(atoms2,x2,box2);
 469
 470   fprintf(stderr,"Select group from second structure\n");
 471   get_index(atoms2,opt2fn_null("-n2",NFILE,fnm),
 472             1,&isize2,&index2,&groupnames2);
 473
 474   if (bName) {
 475     find_matching_names(&isize1, index1, atoms1, &isize2, index2, atoms2);
 476     if (matchndxfile) {
 477       fp = ffopen(matchndxfile,"w");
 478       fprintf(fp, "; Matching atoms between %s from %s and %s from %s\n",
 479               groupnames1, conf1file, groupnames2, conf2file);
 480       fprintf(fp, "[ Match_%s_%s ]\n", conf1file, groupnames1);
 481       for(i=0; i<isize1; i++)
 482         fprintf(fp, "%4u%s", index1[i]+1, (i%15==14 || i==isize1-1)?"\n":" ");
 483       fprintf(fp, "[ Match_%s_%s ]\n", conf2file, groupnames2);
 484       for(i=0; i<isize2; i++)
 485         fprintf(fp, "%4u%s", index2[i]+1, (i%15==14 || i==isize2-1)?"\n":" ");
 486     }
 487   }
 488
 489   /* check isizes, must be equal */
 490   if ( isize2 != isize1 )
 491     gmx_fatal(FARGS,"You selected groups with differen number of atoms.\n");
 492
 493   for(i=0; i<isize1; i++) {
 494     name1=*atoms1->atomname[index1[i]];
 495     name2=*atoms2->atomname[index2[i]];
 496     if (strcmp(name1,name2)) {
 497       if (warn < 20)
 498         fprintf(stderr,
 499                 "Warning: atomnames at index %d don't match: %u %s, %u %s\n",
 500                 i+1,index1[i]+1,name1,index2[i]+1,name2);
 501       warn++;
 502     }
 503     if (!bMW) {
 504       atoms1->atom[index1[i]].m = 1;
 505       atoms2->atom[index2[i]].m = 1;
 506     }
 507   }
 508   if (warn)
 509     fprintf(stderr,"%d atomname%s did not match\n",warn,(warn==1) ? "":"s");
 510
 511   if (bFit) {
 512     /* calculate and remove center of mass of structures */
 513     calc_rm_cm(isize1, index1, atoms1, x1, xcm1);
 514     calc_rm_cm(isize2, index2, atoms2, x2, xcm2);
 515
 516     snew(w_rls,atoms2->nr);
 517     snew(fit_x,atoms2->nr);
 518     for(at=0; (at<isize1); at++) {
 519       w_rls[index2[at]] = atoms1->atom[index1[at]].m;
 520       copy_rvec(x1[index1[at]],fit_x[index2[at]]);
 521     }
 522
 523     /* do the least squares fit to the reference structure */
 524     do_fit(atoms2->nr,w_rls,fit_x,x2);
 525
 526     sfree(fit_x);
 527     sfree(w_rls);
 528     w_rls = NULL;
 529   }
 530   else {
 531     clear_rvec(xcm1);
 532     clear_rvec(xcm2);
 533     w_rls = NULL;
 534   }
 535
 536   /* calculate the rms deviation */
 537   rms     = 0;
 538   totmass = 0;
 539   maxmsd  = -1e18;
 540   minmsd  =  1e18;
 541   snew(msds, isize1);
 542   for(at=0; at<isize1; at++) {
 543     mass = atoms1->atom[index1[at]].m;
 544     for(m=0; m<DIM; m++) {
 545       msd = sqr(x1[index1[at]][m] - x2[index2[at]][m]);
 546       rms += msd*mass;
 547       msds[at] += msd;
 548     }
 549     maxmsd = max(maxmsd, msds[at]);
 550     minmsd = min(minmsd, msds[at]);
 551     totmass += mass;
 552   }
 553   rms = sqrt(rms/totmass);
 554
 555   printf("Root mean square deviation after lsq fit = %g nm\n",rms);
 556   if (bBfac)
 557     printf("Atomic MSD's range from %g to %g nm^2\n",minmsd, maxmsd);
 558
 559   if (bFit) {
 560     /* reset coordinates of reference and fitted structure */
 561     for(i=0; i<atoms1->nr; i++)
 562       for(m=0; m<DIM; m++)
 563         x1[i][m]+=xcm1[m];
 564     for(i=0; i<atoms2->nr; i++)
 565       for(m=0; m<DIM; m++)
 566         x2[i][m]+=xcm1[m];
 567   }
 568
 569   outfile = ftp2fn(efSTO,NFILE,fnm);
 570   switch (fn2ftp(outfile)) {
 571   case efPDB:
 572   case efBRK:
 573   case efENT:
 574     if (bBfac || bLabel) {
 575       srenew(atoms1->pdbinfo, atoms1->nr);
 576       srenew(atoms1->atom, atoms1->nr);      /* Why renew atom? */
 577
 578       /* Avoid segfaults when writing the pdb-file */
 579       for (i=0; i<atoms1->nr; i++) {
 580         atoms1->pdbinfo[i].type = eptAtom;
 581         atoms1->pdbinfo[i].occup = 1.00;
 582         atoms1->pdbinfo[i].bAnisotropic = FALSE;
 583         if (bBfac)
 584           atoms1->pdbinfo[i].bfac = 0;
 585         if (bLabel)
 586           atoms1->resinfo[atoms1->atom[i].resind].chainid = 'A';
 587       }
 588
 589       for(i=0; i<isize1; i++) {
 590         /* atoms1->pdbinfo[index1[i]].type = eptAtom; */
 591 /*      atoms1->pdbinfo[index1[i]].bAnisotropic = FALSE; */
 592         if (bBfac)
 593           atoms1->pdbinfo[index1[i]].bfac = (800*M_PI*M_PI/3.0)*msds[i];
 594 /*      if (bLabel) */
 595 /*        atoms1->resinfo[atoms1->atom[index1[i]].resind].chain = 'A'; */
 596       }
 597       srenew(atoms2->pdbinfo, atoms2->nr);
 598       srenew(atoms2->atom, atoms2->nr);      /* Why renew atom? */
 599
 600       for (i=0; i<atoms2->nr; i++) {
 601         atoms2->pdbinfo[i].type = eptAtom;
 602         atoms2->pdbinfo[i].occup = 1.00;
 603         atoms2->pdbinfo[i].bAnisotropic = FALSE;
 604         if (bBfac)
 605           atoms2->pdbinfo[i].bfac = 0;
 606         if (bLabel)
 607           atoms2->resinfo[atoms1->atom[i].resind].chainid = 'B';
 608       }
 609
 610       for(i=0; i<isize2; i++) {
 611         /* atoms2->pdbinfo[index2[i]].type = eptAtom; */
 612 /*      atoms2->pdbinfo[index2[i]].bAnisotropic = FALSE; */
 613         if (bBfac)
 614           atoms2->pdbinfo[index2[i]].bfac = (800*M_PI*M_PI/3.0)*msds[i];
 615 /*      if (bLabel) */
 616 /*        atoms2->resinfo[atoms2->atom[index2[i]].resind].chain = 'B'; */
 617       }
 618     }
 619     fp=ffopen(outfile,"w");
 620     if (!bOne)
 621       write_pdbfile(fp,title1,atoms1,x1,ePBC1,box1,' ',1,NULL,TRUE);
 622     write_pdbfile(fp,title2,atoms2,x2,ePBC2,box2,' ',bOne ? -1 : 2,NULL,TRUE);
 623     ffclose(fp);
 624     break;
 625   case efGRO:
 626     if (bBfac)
 627       fprintf(stderr,"WARNING: cannot write B-factor values to gro file\n");
 628     fp=ffopen(outfile,"w");
 629     if (!bOne)
 630       write_hconf_p(fp,title1,atoms1,3,x1,v1,box1);
 631     write_hconf_p(fp,title2,atoms2,3,x2,v2,box2);
 632     ffclose(fp);
 633     break;
 634   default:
 635     if (bBfac)
 636       fprintf(stderr,"WARNING: cannot write B-factor values to %s file\n",
 637               ftp2ext(fn2ftp(outfile)));
 638     if (!bOne)
 639       fprintf(stderr,
 640               "WARNING: cannot write the reference structure to %s file\n",
 641               ftp2ext(fn2ftp(outfile)));
 642     write_sto_conf(outfile,title2,atoms2,x2,v2,ePBC2,box2);
 643     break;
 644   }
 645
 646   view_all(oenv,NFILE, fnm);
 647
 648   thanx(stderr);
 649
 650   return 0;
 651 }