src/gromacs/gmxana/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 "gromacs/fileio/filenm.h"
  40 #include "smalloc.h"
  41 #include "macros.h"
  42 #include <math.h>
  43 #include "typedefs.h"
  44 #include "xvgr.h"
  45 #include "statutil.h"
  46 #include "gromacs/fileio/tpxio.h"
  47 #include "string2.h"
  48 #include "vec.h"
  49 #include "index.h"
  50 #include "pbc.h"
  51 #include "gmx_fatal.h"
  52 #include "gromacs/fileio/futil.h"
  53 #include "gromacs/fileio/confio.h"
  54 #include "gromacs/fileio/pdbio.h"
  55 #include "txtdump.h"
  56 #include "do_fit.h"
  57 #include "viewit.h"
  58 #include "rmpbc.h"
  59 #include "gmx_ana.h"
  60
  61
  62 void calc_rm_cm(int isize, atom_id index[], t_atoms *atoms, rvec x[], rvec xcm)
  63 {
  64     int  i, d;
  65     real tm, m;
  66
  67     /* calculate and remove center of mass of reference structure */
  68     tm = 0;
  69     clear_rvec(xcm);
  70     for (i = 0; i < isize; i++)
  71     {
  72         m = atoms->atom[index[i]].m;
  73         for (d = 0; d < DIM; d++)
  74         {
  75             xcm[d] += m*x[index[i]][d];
  76         }
  77         tm += m;
  78     }
  79     svmul(1/tm, xcm, xcm);
  80     for (i = 0; i < atoms->nr; i++)
  81     {
  82         rvec_dec(x[i], xcm);
  83     }
  84 }
  85
  86 int build_res_index(int isize, atom_id index[], t_atom atom[], int rindex[])
  87 {
  88     int i, r;
  89
  90     r         = 0;
  91     rindex[r] = atom[index[0]].resind;
  92     r++;
  93     for (i = 1; i < isize; i++)
  94     {
  95         if (atom[index[i]].resind != rindex[r-1])
  96         {
  97             rindex[r] = atom[index[i]].resind;
  98             r++;
  99         }
 100     }
 101
 102     return r;
 103 }
 104
 105 int find_res_end(int i, int isize, atom_id index[], t_atoms *atoms)
 106 {
 107     int rnr;
 108
 109     rnr = atoms->atom[index[i]].resind;
 110     while (i < isize && atoms->atom[index[i]].resind == rnr)
 111     {
 112         i++;
 113     }
 114     return i;
 115 }
 116
 117 int debug_strcmp(char s1[], char s2[])
 118 {
 119     if (debug)
 120     {
 121         fprintf(debug, " %s-%s", s1, s2);
 122     }
 123     return strcmp(s1, s2);
 124 }
 125
 126 int find_next_match_atoms_in_res(int *i1, atom_id index1[],
 127                                  int m1, char **atnms1[],
 128                                  int *i2, atom_id index2[],
 129                                  int m2, char **atnms2[])
 130 {
 131     int      dx, dy, dmax, cmp;
 132     gmx_bool bFW = FALSE;
 133
 134     dx   = dy = 0;
 135     cmp  = NOTSET;
 136     dmax = max(m1-*i1, m2-*i2);
 137     for (dx = 0; dx < dmax && cmp != 0; dx++)
 138     {
 139         for (dy = dx; dy < dmax && cmp != 0; dy++)
 140         {
 141             if (dx || dy)
 142             {
 143                 if (debug)
 144                 {
 145                     fprintf(debug, ".");
 146                 }
 147                 cmp = NOTSET;
 148                 if (*i1+dx < m1 && *i2+dy < m2)
 149                 {
 150                     bFW = TRUE;
 151                     cmp = debug_strcmp(*atnms1[index1[*i1+dx]], *atnms2[index2[*i2+dy]]);
 152                     if (debug)
 153                     {
 154                         fprintf(debug, "(%d %d)", *i1+dx, *i2+dy);
 155                     }
 156                 }
 157                 if (cmp != 0 && *i1+dy < m1 && *i2+dx < m2)
 158                 {
 159                     bFW = FALSE;
 160                     cmp = debug_strcmp(*atnms1[index1[*i1+dy]], *atnms2[index2[*i2+dx]]);
 161                     if (debug)
 162                     {
 163                         fprintf(debug, "(%d %d)", *i1+dy, *i2+dx);
 164                     }
 165                 }
 166             }
 167         }
 168     }
 169     /* apparently, dx and dy are incremented one more time
 170        as the loops terminate: we correct this here */
 171     dx--;
 172     dy--;
 173     if (cmp == 0)
 174     {
 175         if (debug)
 176         {
 177             fprintf(debug, "{%d %d}", *i1+bFW ? dx : dy, *i2+bFW ? dy : dx );
 178         }
 179         if (bFW)
 180         {
 181             *i1 += dx;
 182             *i2 += dy;
 183         }
 184         else
 185         {
 186             *i1 += dy;
 187             *i2 += dx;
 188         }
 189     }
 190
 191     return cmp;
 192 }
 193
 194 static int find_next_match_res(int *rnr1, int isize1,
 195                                int index1[], t_resinfo *resinfo1,
 196                                int *rnr2, int isize2,
 197                                int index2[], t_resinfo *resinfo2)
 198 {
 199     int      dx, dy, dmax, cmp, rr1, rr2;
 200     gmx_bool bFW = FALSE, bFF = FALSE;
 201
 202     dx  = dy = 0;
 203     rr1 = 0;
 204     while (rr1 < isize1 && *rnr1 != index1[rr1])
 205     {
 206         rr1++;
 207     }
 208     rr2 = 0;
 209     while (rr2 < isize2 && *rnr2 != index2[rr2])
 210     {
 211         rr2++;
 212     }
 213
 214     cmp  = NOTSET;
 215     dmax = max(isize1-rr1, isize2-rr2);
 216     if (debug)
 217     {
 218         fprintf(debug, " R:%d-%d:%d-%d:%d ",
 219                 rr1, isize1, rr2, isize2, dmax);
 220     }
 221     for (dx = 0; dx < dmax && cmp != 0; dx++)
 222     {
 223         for (dy = 0; dy <= dx && cmp != 0; dy++)
 224         {
 225             if (dx != dy)
 226             {
 227                 cmp = NOTSET;
 228                 if (rr1+dx < isize1 && rr2+dy < isize2)
 229                 {
 230                     bFW = TRUE;
 231                     cmp = debug_strcmp(*resinfo1[index1[rr1+dx]].name,
 232                                        *resinfo2[index2[rr2+dy]].name);
 233                     if (debug)
 234                     {
 235                         fprintf(debug, "(%d %d)", rr1+dx, rr2+dy);
 236                     }
 237                 }
 238                 if (cmp != 0 && rr1+dy < isize1 && rr2+dx < isize2)
 239                 {
 240                     bFW = FALSE;
 241                     cmp = debug_strcmp(*resinfo1[index1[rr1+dy]].name,
 242                                        *resinfo2[index2[rr2+dx]].name);
 243                     if (debug)
 244                     {
 245                         fprintf(debug, "(%d %d)", rr1+dy, rr2+dx);
 246                     }
 247                 }
 248                 if (dx != 0 && cmp != 0 && rr1+dx < isize1 && rr2+dx < isize2)
 249                 {
 250                     bFF = TRUE;
 251                     cmp = debug_strcmp(*resinfo1[index1[rr1+dx]].name,
 252                                        *resinfo2[index2[rr2+dx]].name);
 253                     if (debug)
 254                     {
 255                         fprintf(debug, "(%d %d)", rr1+dx, rr2+dx);
 256                     }
 257                 }
 258                 else
 259                 {
 260                     bFF = FALSE;
 261                 }
 262             }
 263         }
 264     }
 265     /* apparently, dx and dy are incremented one more time
 266        as the loops terminate: we correct this here */
 267     dx--;
 268     dy--;
 269     /* if we skipped equal on both sides, only skip one residue
 270        to allow for single mutations of residues... */
 271     if (bFF)
 272     {
 273         if (debug)
 274         {
 275             fprintf(debug, "%d.%d.%dX%sX%s", dx, rr1, rr2,
 276                     *resinfo1[index1[rr1+1]].name,
 277                     *resinfo2[index2[rr2+1]].name);
 278         }
 279         dx = 1;
 280     }
 281     if (cmp == 0)
 282     {
 283         if (debug)
 284         {
 285             fprintf(debug, "!");
 286         }
 287         if (bFF)
 288         {
 289             rr1 += dx;
 290             rr2 += dx;
 291         }
 292         else
 293         if (bFW)
 294         {
 295             rr1 += dx;
 296             rr2 += dy;
 297         }
 298         else
 299         {
 300             rr1 += dy;
 301             rr2 += dx;
 302         }
 303         *rnr1 = index1[rr1];
 304         *rnr2 = index2[rr2];
 305     }
 306
 307     return cmp;
 308 }
 309
 310 int find_first_atom_in_res(int rnr, int isize, atom_id index[], t_atom atom[])
 311 {
 312     int i;
 313
 314     i = 0;
 315     while (i < isize && atom[index[i]].resind != rnr)
 316     {
 317         i++;
 318     }
 319
 320     if (atom[index[i]].resind == rnr)
 321     {
 322         return i;
 323     }
 324     else
 325     {
 326         return NOTSET;
 327     }
 328 }
 329
 330 void find_matching_names(int *isize1, atom_id index1[], t_atoms *atoms1,
 331                          int *isize2, atom_id index2[], t_atoms *atoms2)
 332 {
 333     int        i, i1, i2, ii1, ii2, m1, m2;
 334     int        atcmp, rescmp;
 335     int        r, rnr1, rnr2, prnr1, prnr2;
 336     int        rsize1, rsize2;
 337     int       *rindex1, *rindex2;
 338     char      *resnm1, *resnm2, *atnm1, *atnm2;
 339     char    ***atnms1, ***atnms2;
 340     t_resinfo *resinfo1, *resinfo2;
 341
 342     /* set some handy shortcuts */
 343     resinfo1 = atoms1->resinfo;
 344     atnms1   = atoms1->atomname;
 345     resinfo2 = atoms2->resinfo;
 346     atnms2   = atoms2->atomname;
 347
 348     /* build indexes of selected residues */
 349     snew(rindex1, atoms1->nres);
 350     rsize1 = build_res_index(*isize1, index1, atoms1->atom, rindex1);
 351     snew(rindex2, atoms2->nres);
 352     rsize2 = build_res_index(*isize2, index2, atoms2->atom, rindex2);
 353
 354     i1    = i2 = 0;
 355     ii1   = ii2 = 0;
 356     atcmp = rescmp = 0;
 357     prnr1 = prnr2 = NOTSET;
 358     if (debug)
 359     {
 360         fprintf(debug, "Find matching names: %d, %d\n", *isize1, *isize2);
 361     }
 362     while (atcmp == 0 && i1 < *isize1 && i2 < *isize2)
 363     {
 364         /* prologue */
 365         rnr1 = atoms1->atom[index1[i1]].resind;
 366         rnr2 = atoms2->atom[index2[i2]].resind;
 367         if (rnr1 != prnr1 || rnr2 != prnr2)
 368         {
 369             if (debug)
 370             {
 371                 fprintf(debug, "R: %s%d %s%d\n",
 372                         *resinfo1[rnr1].name, rnr1, *resinfo2[rnr2].name, rnr2);
 373             }
 374             rescmp = strcmp(*resinfo1[rnr1].name, *resinfo2[rnr2].name);
 375         }
 376         if (debug)
 377         {
 378             fprintf(debug, "comparing %d %d", i1, i2);
 379         }
 380         atcmp = debug_strcmp(*atnms1[index1[i1]], *atnms2[index2[i2]]);
 381
 382         /* the works */
 383         if (atcmp != 0) /* no match -> find match within residues */
 384         {
 385             m1 = find_res_end(i1, *isize1, index1, atoms1);
 386             m2 = find_res_end(i2, *isize2, index2, atoms2);
 387             if (debug)
 388             {
 389                 fprintf(debug, " [%d<%d %d<%d]", i1, m1, i2, m2);
 390             }
 391             atcmp = find_next_match_atoms_in_res(&i1, index1, m1, atnms1,
 392                                                  &i2, index2, m2, atnms2);
 393             if (debug)
 394             {
 395                 fprintf(debug, " -> %d %d %s-%s", i1, i2,
 396                         *atnms1[index1[i1]], *atnms2[index2[i2]]);
 397             }
 398
 399         }
 400         if (atcmp != 0) /* still no match -> next residue(s) */
 401         {
 402             prnr1  = rnr1;
 403             prnr2  = rnr2;
 404             rescmp = find_next_match_res(&rnr1, rsize1, rindex1, resinfo1,
 405                                          &rnr2, rsize2, rindex2, resinfo2);
 406             if (rnr1 != prnr1)
 407             {
 408                 i1 = find_first_atom_in_res(rnr1, *isize1, index1, atoms1->atom);
 409             }
 410             if (rnr2 != prnr2)
 411             {
 412                 i2 = find_first_atom_in_res(rnr2, *isize2, index2, atoms2->atom);
 413             }
 414             if (debug)
 415             {
 416                 fprintf(debug, " -> %s%d-%s%d %s%d-%s%d",
 417                         *resinfo1[rnr1].name, rnr1,
 418                         *atnms1[index1[i1]], index1[i1],
 419                         *resinfo2[rnr2].name, rnr2,
 420                         *atnms2[index2[i2]], index2[i2]);
 421             }
 422             m1 = find_res_end(i1, *isize1, index1, atoms1);
 423             m2 = find_res_end(i2, *isize2, index2, atoms2);
 424             if (debug)
 425             {
 426                 fprintf(debug, " [%d<%d %d<%d]", i1, m1, i2, m2);
 427             }
 428             atcmp = find_next_match_atoms_in_res(&i1, index1, m1, atnms1,
 429                                                  &i2, index2, m2, atnms2);
 430             if (debug)
 431             {
 432                 fprintf(debug, " -> %d %d %s-%s", i1, i2,
 433                         *atnms1[index1[i1]], *atnms2[index2[i2]]);
 434             }
 435         }
 436         if (debug)
 437         {
 438             fprintf(debug, "(%d %d): %d %d\n", ii1, ii2, atcmp, rescmp);
 439         }
 440         if (atcmp == 0) /* if match -> store indices */
 441         {
 442             index1[ii1++] = index1[i1];
 443             index2[ii2++] = index2[i2];
 444         }
 445         i1++;
 446         i2++;
 447
 448         /* epilogue */
 449         prnr1 = rnr1;
 450         prnr2 = rnr2;
 451     }
 452
 453     if (ii1 != ii2)
 454     {
 455         gmx_fatal(FARGS, "DEATH HORROR: non-equal number of matching atoms!\n");
 456     }
 457     if (ii1 == i1 && ii2 == i2)
 458     {
 459         printf("All atoms in index groups 1 and 2 match\n");
 460     }
 461     else
 462     {
 463         if (i1 == i2 && ii1 == ii2)
 464         {
 465             printf("Both index groups modified from %d to %d atoms\n", i1, ii1);
 466         }
 467         else
 468         {
 469             if (ii1 != i1)
 470             {
 471                 printf("Index group 1 modified from %d to %d atoms\n", i1, ii1);
 472             }
 473             if (ii2 != i2)
 474             {
 475                 printf("Index group 2 modified from %d to %d atoms\n", i2, ii2);
 476             }
 477         }
 478         *isize1 = ii1;
 479         *isize2 = ii2;
 480     }
 481 }
 482 /* 1 */
 483
 484 int gmx_confrms(int argc, char *argv[])
 485 {
 486     const char     *desc[] = {
 487         "[TT]g_confrms[tt] computes the root mean square deviation (RMSD) of two",
 488         "structures after least-squares fitting the second structure on the first one.",
 489         "The two structures do NOT need to have the same number of atoms,",
 490         "only the two index groups used for the fit need to be identical.",
 491         "With [TT]-name[tt] only matching atom names from the selected groups",
 492         "will be used for the fit and RMSD calculation. This can be useful ",
 493         "when comparing mutants of a protein.",
 494         "[PAR]",
 495         "The superimposed structures are written to file. In a [TT].pdb[tt] file",
 496         "the two structures will be written as separate models",
 497         "(use [TT]rasmol -nmrpdb[tt]). Also in a [TT].pdb[tt] file, B-factors",
 498         "calculated from the atomic MSD values can be written with [TT]-bfac[tt].",
 499     };
 500     static gmx_bool bOne  = FALSE, bRmpbc = FALSE, bMW = TRUE, bName = FALSE,
 501                     bBfac = FALSE, bFit = TRUE, bLabel = FALSE;
 502
 503     t_pargs  pa[] = {
 504         { "-one", FALSE, etBOOL, {&bOne},   "Only write the fitted structure to file" },
 505         { "-mw",  FALSE, etBOOL, {&bMW},    "Mass-weighted fitting and RMSD" },
 506         { "-pbc", FALSE, etBOOL, {&bRmpbc}, "Try to make molecules whole again" },
 507         { "-fit", FALSE, etBOOL, {&bFit},
 508           "Do least squares superposition of the target structure to the reference" },
 509         { "-name", FALSE, etBOOL, {&bName},
 510           "Only compare matching atom names" },
 511         { "-label", FALSE, etBOOL, {&bLabel},
 512           "Added chain labels A for first and B for second structure"},
 513         { "-bfac", FALSE, etBOOL, {&bBfac},
 514           "Output B-factors from atomic MSD values" }
 515     };
 516     t_filenm fnm[] = {
 517         { efTPS, "-f1",  "conf1.gro", ffREAD  },
 518         { efSTX, "-f2",  "conf2",     ffREAD  },
 519         { efSTO, "-o",   "fit.pdb",   ffWRITE },
 520         { efNDX, "-n1", "fit1.ndx",  ffOPTRD },
 521         { efNDX, "-n2", "fit2.ndx",  ffOPTRD },
 522         { efNDX, "-no", "match.ndx", ffOPTWR }
 523     };
 524 #define NFILE asize(fnm)
 525
 526     /* the two structure files */
 527     const char  *conf1file, *conf2file, *matchndxfile, *outfile;
 528     FILE        *fp;
 529     char         title1[STRLEN], title2[STRLEN], *name1, *name2;
 530     t_topology  *top1, *top2;
 531     int          ePBC1, ePBC2;
 532     t_atoms     *atoms1, *atoms2;
 533     int          warn = 0;
 534     atom_id      at;
 535     real        *w_rls, mass, totmass;
 536     rvec        *x1, *v1, *x2, *v2, *fit_x;
 537     matrix       box1, box2;
 538
 539     output_env_t oenv;
 540
 541     /* counters */
 542     int     i, j, m;
 543
 544     /* center of mass calculation */
 545     real    tmas1, tmas2;
 546     rvec    xcm1, xcm2;
 547
 548     /* variables for fit */
 549     char    *groupnames1, *groupnames2;
 550     int      isize1, isize2;
 551     atom_id *index1, *index2;
 552     real     rms, msd, minmsd, maxmsd;
 553     real    *msds;
 554
 555
 556     if (!parse_common_args(&argc, argv, PCA_BE_NICE | PCA_CAN_VIEW,
 557                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 558     {
 559         return 0;
 560     }
 561     matchndxfile = opt2fn_null("-no", NFILE, fnm);
 562     conf1file    = ftp2fn(efTPS, NFILE, fnm);
 563     conf2file    = ftp2fn(efSTX, NFILE, fnm);
 564
 565     /* reading reference structure from first structure file */
 566     fprintf(stderr, "\nReading first structure file\n");
 567     snew(top1, 1);
 568     read_tps_conf(conf1file, title1, top1, &ePBC1, &x1, &v1, box1, TRUE);
 569     atoms1 = &(top1->atoms);
 570     fprintf(stderr, "%s\nContaining %d atoms in %d residues\n",
 571             title1, atoms1->nr, atoms1->nres);
 572
 573     if (bRmpbc)
 574     {
 575         rm_gropbc(atoms1, x1, box1);
 576     }
 577
 578     fprintf(stderr, "Select group from first structure\n");
 579     get_index(atoms1, opt2fn_null("-n1", NFILE, fnm),
 580               1, &isize1, &index1, &groupnames1);
 581     printf("\n");
 582
 583     if (bFit && (isize1 < 3))
 584     {
 585         gmx_fatal(FARGS, "Need >= 3 points to fit!\n");
 586     }
 587
 588     /* reading second structure file */
 589     fprintf(stderr, "\nReading second structure file\n");
 590     snew(top2, 1);
 591     read_tps_conf(conf2file, title2, top2, &ePBC2, &x2, &v2, box2, TRUE);
 592     atoms2 = &(top2->atoms);
 593     fprintf(stderr, "%s\nContaining %d atoms in %d residues\n",
 594             title2, atoms2->nr, atoms2->nres);
 595
 596     if (bRmpbc)
 597     {
 598         rm_gropbc(atoms2, x2, box2);
 599     }
 600
 601     fprintf(stderr, "Select group from second structure\n");
 602     get_index(atoms2, opt2fn_null("-n2", NFILE, fnm),
 603               1, &isize2, &index2, &groupnames2);
 604
 605     if (bName)
 606     {
 607         find_matching_names(&isize1, index1, atoms1, &isize2, index2, atoms2);
 608         if (matchndxfile)
 609         {
 610             fp = ffopen(matchndxfile, "w");
 611             fprintf(fp, "; Matching atoms between %s from %s and %s from %s\n",
 612                     groupnames1, conf1file, groupnames2, conf2file);
 613             fprintf(fp, "[ Match_%s_%s ]\n", conf1file, groupnames1);
 614             for (i = 0; i < isize1; i++)
 615             {
 616                 fprintf(fp, "%4u%s", index1[i]+1, (i%15 == 14 || i == isize1-1) ? "\n" : " ");
 617             }
 618             fprintf(fp, "[ Match_%s_%s ]\n", conf2file, groupnames2);
 619             for (i = 0; i < isize2; i++)
 620             {
 621                 fprintf(fp, "%4u%s", index2[i]+1, (i%15 == 14 || i == isize2-1) ? "\n" : " ");
 622             }
 623         }
 624     }
 625
 626     /* check isizes, must be equal */
 627     if (isize2 != isize1)
 628     {
 629         gmx_fatal(FARGS, "You selected groups with differen number of atoms.\n");
 630     }
 631
 632     for (i = 0; i < isize1; i++)
 633     {
 634         name1 = *atoms1->atomname[index1[i]];
 635         name2 = *atoms2->atomname[index2[i]];
 636         if (strcmp(name1, name2))
 637         {
 638             if (warn < 20)
 639             {
 640                 fprintf(stderr,
 641                         "Warning: atomnames at index %d don't match: %u %s, %u %s\n",
 642                         i+1, index1[i]+1, name1, index2[i]+1, name2);
 643             }
 644             warn++;
 645         }
 646         if (!bMW)
 647         {
 648             atoms1->atom[index1[i]].m = 1;
 649             atoms2->atom[index2[i]].m = 1;
 650         }
 651     }
 652     if (warn)
 653     {
 654         fprintf(stderr, "%d atomname%s did not match\n", warn, (warn == 1) ? "" : "s");
 655     }
 656
 657     if (bFit)
 658     {
 659         /* calculate and remove center of mass of structures */
 660         calc_rm_cm(isize1, index1, atoms1, x1, xcm1);
 661         calc_rm_cm(isize2, index2, atoms2, x2, xcm2);
 662
 663         snew(w_rls, atoms2->nr);
 664         snew(fit_x, atoms2->nr);
 665         for (at = 0; (at < isize1); at++)
 666         {
 667             w_rls[index2[at]] = atoms1->atom[index1[at]].m;
 668             copy_rvec(x1[index1[at]], fit_x[index2[at]]);
 669         }
 670
 671         /* do the least squares fit to the reference structure */
 672         do_fit(atoms2->nr, w_rls, fit_x, x2);
 673
 674         sfree(fit_x);
 675         sfree(w_rls);
 676         w_rls = NULL;
 677     }
 678     else
 679     {
 680         clear_rvec(xcm1);
 681         clear_rvec(xcm2);
 682         w_rls = NULL;
 683     }
 684
 685     /* calculate the rms deviation */
 686     rms     = 0;
 687     totmass = 0;
 688     maxmsd  = -1e18;
 689     minmsd  =  1e18;
 690     snew(msds, isize1);
 691     for (at = 0; at < isize1; at++)
 692     {
 693         mass = atoms1->atom[index1[at]].m;
 694         for (m = 0; m < DIM; m++)
 695         {
 696             msd       = sqr(x1[index1[at]][m] - x2[index2[at]][m]);
 697             rms      += msd*mass;
 698             msds[at] += msd;
 699         }
 700         maxmsd   = max(maxmsd, msds[at]);
 701         minmsd   = min(minmsd, msds[at]);
 702         totmass += mass;
 703     }
 704     rms = sqrt(rms/totmass);
 705
 706     printf("Root mean square deviation after lsq fit = %g nm\n", rms);
 707     if (bBfac)
 708     {
 709         printf("Atomic MSD's range from %g to %g nm^2\n", minmsd, maxmsd);
 710     }
 711
 712     if (bFit)
 713     {
 714         /* reset coordinates of reference and fitted structure */
 715         for (i = 0; i < atoms1->nr; i++)
 716         {
 717             for (m = 0; m < DIM; m++)
 718             {
 719                 x1[i][m] += xcm1[m];
 720             }
 721         }
 722         for (i = 0; i < atoms2->nr; i++)
 723         {
 724             for (m = 0; m < DIM; m++)
 725             {
 726                 x2[i][m] += xcm1[m];
 727             }
 728         }
 729     }
 730
 731     outfile = ftp2fn(efSTO, NFILE, fnm);
 732     switch (fn2ftp(outfile))
 733     {
 734         case efPDB:
 735         case efBRK:
 736         case efENT:
 737             if (bBfac || bLabel)
 738             {
 739                 srenew(atoms1->pdbinfo, atoms1->nr);
 740                 srenew(atoms1->atom, atoms1->nr); /* Why renew atom? */
 741
 742                 /* Avoid segfaults when writing the pdb-file */
 743                 for (i = 0; i < atoms1->nr; i++)
 744                 {
 745                     atoms1->pdbinfo[i].type         = eptAtom;
 746                     atoms1->pdbinfo[i].occup        = 1.00;
 747                     atoms1->pdbinfo[i].bAnisotropic = FALSE;
 748                     if (bBfac)
 749                     {
 750                         atoms1->pdbinfo[i].bfac = 0;
 751                     }
 752                     if (bLabel)
 753                     {
 754                         atoms1->resinfo[atoms1->atom[i].resind].chainid = 'A';
 755                     }
 756                 }
 757
 758                 for (i = 0; i < isize1; i++)
 759                 {
 760                     /* atoms1->pdbinfo[index1[i]].type = eptAtom; */
 761 /*  atoms1->pdbinfo[index1[i]].bAnisotropic = FALSE; */
 762                     if (bBfac)
 763                     {
 764                         atoms1->pdbinfo[index1[i]].bfac = (800*M_PI*M_PI/3.0)*msds[i];
 765                     }
 766 /*  if (bLabel) */
 767 /*    atoms1->resinfo[atoms1->atom[index1[i]].resind].chain = 'A'; */
 768                 }
 769                 srenew(atoms2->pdbinfo, atoms2->nr);
 770                 srenew(atoms2->atom, atoms2->nr); /* Why renew atom? */
 771
 772                 for (i = 0; i < atoms2->nr; i++)
 773                 {
 774                     atoms2->pdbinfo[i].type         = eptAtom;
 775                     atoms2->pdbinfo[i].occup        = 1.00;
 776                     atoms2->pdbinfo[i].bAnisotropic = FALSE;
 777                     if (bBfac)
 778                     {
 779                         atoms2->pdbinfo[i].bfac = 0;
 780                     }
 781                     if (bLabel)
 782                     {
 783                         atoms2->resinfo[atoms1->atom[i].resind].chainid = 'B';
 784                     }
 785                 }
 786
 787                 for (i = 0; i < isize2; i++)
 788                 {
 789                     /* atoms2->pdbinfo[index2[i]].type = eptAtom; */
 790 /*  atoms2->pdbinfo[index2[i]].bAnisotropic = FALSE; */
 791                     if (bBfac)
 792                     {
 793                         atoms2->pdbinfo[index2[i]].bfac = (800*M_PI*M_PI/3.0)*msds[i];
 794                     }
 795 /*  if (bLabel) */
 796 /*    atoms2->resinfo[atoms2->atom[index2[i]].resind].chain = 'B'; */
 797                 }
 798             }
 799             fp = ffopen(outfile, "w");
 800             if (!bOne)
 801             {
 802                 write_pdbfile(fp, title1, atoms1, x1, ePBC1, box1, ' ', 1, NULL, TRUE);
 803             }
 804             write_pdbfile(fp, title2, atoms2, x2, ePBC2, box2, ' ', bOne ? -1 : 2, NULL, TRUE);
 805             ffclose(fp);
 806             break;
 807         case efGRO:
 808             if (bBfac)
 809             {
 810                 fprintf(stderr, "WARNING: cannot write B-factor values to gro file\n");
 811             }
 812             fp = ffopen(outfile, "w");
 813             if (!bOne)
 814             {
 815                 write_hconf_p(fp, title1, atoms1, 3, x1, v1, box1);
 816             }
 817             write_hconf_p(fp, title2, atoms2, 3, x2, v2, box2);
 818             ffclose(fp);
 819             break;
 820         default:
 821             if (bBfac)
 822             {
 823                 fprintf(stderr, "WARNING: cannot write B-factor values to %s file\n",
 824                         ftp2ext(fn2ftp(outfile)));
 825             }
 826             if (!bOne)
 827             {
 828                 fprintf(stderr,
 829                         "WARNING: cannot write the reference structure to %s file\n",
 830                         ftp2ext(fn2ftp(outfile)));
 831             }
 832             write_sto_conf(outfile, title2, atoms2, x2, v2, ePBC2, box2);
 833             break;
 834     }
 835
 836     view_all(oenv, NFILE, fnm);
 837
 838     return 0;
 839 }