src/tools/gmx_confrms.c

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