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