src/gromacs/gmxana/gmx_confrms.cpp

   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,2015,2016,2017, 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 <cmath>
  40 #include <cstring>
  41
  42 #include <algorithm>
  43
  44 #include "gromacs/commandline/filenm.h"
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/commandline/viewit.h"
  47 #include "gromacs/fileio/confio.h"
  48 #include "gromacs/fileio/groio.h"
  49 #include "gromacs/fileio/pdbio.h"
  50 #include "gromacs/gmxana/gmx_ana.h"
  51 #include "gromacs/math/do_fit.h"
  52 #include "gromacs/math/functions.h"
  53 #include "gromacs/math/utilities.h"
  54 #include "gromacs/math/vec.h"
  55 #include "gromacs/pbcutil/rmpbc.h"
  56 #include "gromacs/topology/index.h"
  57 #include "gromacs/topology/topology.h"
  58 #include "gromacs/utility/arraysize.h"
  59 #include "gromacs/utility/cstringutil.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/futil.h"
  62 #include "gromacs/utility/smalloc.h"
  63
  64 static const int NOTSET = -9368163;
  65
  66 static void calc_rm_cm(int isize, int index[], const 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     {
  76         m = atoms->atom[index[i]].m;
  77         for (d = 0; d < DIM; d++)
  78         {
  79             xcm[d] += m*x[index[i]][d];
  80         }
  81         tm += m;
  82     }
  83     svmul(1/tm, xcm, xcm);
  84     for (i = 0; i < atoms->nr; i++)
  85     {
  86         rvec_dec(x[i], xcm);
  87     }
  88 }
  89
  90 static int build_res_index(int isize, int index[], t_atom atom[], int rindex[])
  91 {
  92     int i, r;
  93
  94     r         = 0;
  95     rindex[r] = atom[index[0]].resind;
  96     r++;
  97     for (i = 1; i < isize; i++)
  98     {
  99         if (atom[index[i]].resind != rindex[r-1])
 100         {
 101             rindex[r] = atom[index[i]].resind;
 102             r++;
 103         }
 104     }
 105
 106     return r;
 107 }
 108
 109 static int find_res_end(int i, int isize, int index[], const t_atoms *atoms)
 110 {
 111     int rnr;
 112
 113     rnr = atoms->atom[index[i]].resind;
 114     while (i < isize && atoms->atom[index[i]].resind == rnr)
 115     {
 116         i++;
 117     }
 118     return i;
 119 }
 120
 121 static int debug_strcmp(char s1[], char s2[])
 122 {
 123     if (debug)
 124     {
 125         fprintf(debug, " %s-%s", s1, s2);
 126     }
 127     return std::strcmp(s1, s2);
 128 }
 129
 130 static int find_next_match_atoms_in_res(int *i1, int index1[],
 131                                         int m1, char **atnms1[],
 132                                         int *i2, int index2[],
 133                                         int m2, char **atnms2[])
 134 {
 135     int      dx, dy, dmax, cmp;
 136     gmx_bool bFW = FALSE;
 137
 138     cmp  = NOTSET;
 139     dmax = std::max(m1-*i1, m2-*i2);
 140     for (dx = 0, dy = 0; dx < dmax && cmp != 0; dx++)
 141     {
 142         for (dy = dx; dy < dmax && cmp != 0; dy++)
 143         {
 144             if (dx || dy)
 145             {
 146                 if (debug)
 147                 {
 148                     fprintf(debug, ".");
 149                 }
 150                 cmp = NOTSET;
 151                 if (*i1+dx < m1 && *i2+dy < m2)
 152                 {
 153                     bFW = TRUE;
 154                     cmp = debug_strcmp(*atnms1[index1[*i1+dx]], *atnms2[index2[*i2+dy]]);
 155                     if (debug)
 156                     {
 157                         fprintf(debug, "(%d %d)", *i1+dx, *i2+dy);
 158                     }
 159                 }
 160                 if (cmp != 0 && *i1+dy < m1 && *i2+dx < m2)
 161                 {
 162                     bFW = FALSE;
 163                     cmp = debug_strcmp(*atnms1[index1[*i1+dy]], *atnms2[index2[*i2+dx]]);
 164                     if (debug)
 165                     {
 166                         fprintf(debug, "(%d %d)", *i1+dy, *i2+dx);
 167                     }
 168                 }
 169             }
 170         }
 171     }
 172     /* apparently, dx and dy are incremented one more time
 173        as the loops terminate: we correct this here */
 174     dx--;
 175     dy--;
 176     if (cmp == 0)
 177     {
 178         if (debug)
 179         {
 180             fprintf(debug, "{%d %d}", *i1 + (bFW ? dx : dy), *i2 + (bFW ? dy : dx) );
 181         }
 182         if (bFW)
 183         {
 184             *i1 += dx;
 185             *i2 += dy;
 186         }
 187         else
 188         {
 189             *i1 += dy;
 190             *i2 += dx;
 191         }
 192     }
 193
 194     return cmp;
 195 }
 196
 197 static int find_next_match_res(int *rnr1, int isize1,
 198                                int index1[], t_resinfo *resinfo1,
 199                                int *rnr2, int isize2,
 200                                int index2[], t_resinfo *resinfo2)
 201 {
 202     int      dx, dy, dmax, cmp, rr1, rr2;
 203     gmx_bool bFW = FALSE, bFF = FALSE;
 204
 205     rr1 = 0;
 206     while (rr1 < isize1 && *rnr1 != index1[rr1])
 207     {
 208         rr1++;
 209     }
 210     rr2 = 0;
 211     while (rr2 < isize2 && *rnr2 != index2[rr2])
 212     {
 213         rr2++;
 214     }
 215
 216     cmp  = NOTSET;
 217     dmax = std::max(isize1-rr1, isize2-rr2);
 218     if (debug)
 219     {
 220         fprintf(debug, " R:%d-%d:%d-%d:%d ",
 221                 rr1, isize1, rr2, isize2, dmax);
 222     }
 223     for (dx = 0; dx < dmax && cmp != 0; dx++)
 224     {
 225         for (dy = 0; dy <= dx && cmp != 0; dy++)
 226         {
 227             if (dx != dy)
 228             {
 229                 cmp = NOTSET;
 230                 if (rr1+dx < isize1 && rr2+dy < isize2)
 231                 {
 232                     bFW = TRUE;
 233                     cmp = debug_strcmp(*resinfo1[index1[rr1+dx]].name,
 234                                        *resinfo2[index2[rr2+dy]].name);
 235                     if (debug)
 236                     {
 237                         fprintf(debug, "(%d %d)", rr1+dx, rr2+dy);
 238                     }
 239                 }
 240                 if (cmp != 0 && rr1+dy < isize1 && rr2+dx < isize2)
 241                 {
 242                     bFW = FALSE;
 243                     cmp = debug_strcmp(*resinfo1[index1[rr1+dy]].name,
 244                                        *resinfo2[index2[rr2+dx]].name);
 245                     if (debug)
 246                     {
 247                         fprintf(debug, "(%d %d)", rr1+dy, rr2+dx);
 248                     }
 249                 }
 250                 if (dx != 0 && cmp != 0 && rr1+dx < isize1 && rr2+dx < isize2)
 251                 {
 252                     bFF = TRUE;
 253                     cmp = debug_strcmp(*resinfo1[index1[rr1+dx]].name,
 254                                        *resinfo2[index2[rr2+dx]].name);
 255                     if (debug)
 256                     {
 257                         fprintf(debug, "(%d %d)", rr1+dx, rr2+dx);
 258                     }
 259                 }
 260                 else
 261                 {
 262                     bFF = FALSE;
 263                 }
 264             }
 265         }
 266     }
 267     /* apparently, dx and dy are incremented one more time
 268        as the loops terminate: we correct this here */
 269     dx--;
 270     dy--;
 271     /* if we skipped equal on both sides, only skip one residue
 272        to allow for single mutations of residues... */
 273     if (bFF)
 274     {
 275         if (debug)
 276         {
 277             fprintf(debug, "%d.%d.%dX%sX%s", dx, rr1, rr2,
 278                     *resinfo1[index1[rr1+1]].name,
 279                     *resinfo2[index2[rr2+1]].name);
 280         }
 281         dx = 1;
 282     }
 283     if (cmp == 0)
 284     {
 285         if (debug)
 286         {
 287             fprintf(debug, "!");
 288         }
 289         if (bFF)
 290         {
 291             rr1 += dx;
 292             rr2 += dx;
 293         }
 294         else
 295         if (bFW)
 296         {
 297             rr1 += dx;
 298             rr2 += dy;
 299         }
 300         else
 301         {
 302             rr1 += dy;
 303             rr2 += dx;
 304         }
 305         *rnr1 = index1[rr1];
 306         *rnr2 = index2[rr2];
 307     }
 308
 309     return cmp;
 310 }
 311
 312 static int find_first_atom_in_res(int rnr, int isize, int index[], t_atom atom[])
 313 {
 314     int i;
 315
 316     i = 0;
 317     while (i < isize && atom[index[i]].resind != rnr)
 318     {
 319         i++;
 320     }
 321
 322     if (atom[index[i]].resind == rnr)
 323     {
 324         return i;
 325     }
 326     else
 327     {
 328         return NOTSET;
 329     }
 330 }
 331
 332 static void find_matching_names(int *isize1, int index1[], const t_atoms *atoms1,
 333                                 int *isize2, int index2[], const t_atoms *atoms2)
 334 {
 335     int        i1, i2, ii1, ii2, m1, m2;
 336     int        atcmp, rescmp;
 337     int        rnr1, rnr2, prnr1, prnr2;
 338     int        rsize1, rsize2;
 339     int       *rindex1, *rindex2;
 340     char    ***atnms1, ***atnms2;
 341     t_resinfo *resinfo1, *resinfo2;
 342
 343     /* set some handy shortcuts */
 344     resinfo1 = atoms1->resinfo;
 345     atnms1   = atoms1->atomname;
 346     resinfo2 = atoms2->resinfo;
 347     atnms2   = atoms2->atomname;
 348
 349     /* build indexes of selected residues */
 350     snew(rindex1, atoms1->nres);
 351     rsize1 = build_res_index(*isize1, index1, atoms1->atom, rindex1);
 352     snew(rindex2, atoms2->nres);
 353     rsize2 = build_res_index(*isize2, index2, atoms2->atom, rindex2);
 354
 355     i1    = i2 = 0;
 356     ii1   = ii2 = 0;
 357     atcmp = rescmp = 0;
 358     prnr1 = prnr2 = NOTSET;
 359     if (debug)
 360     {
 361         fprintf(debug, "Find matching names: %d, %d\n", *isize1, *isize2);
 362     }
 363     while (atcmp == 0 && i1 < *isize1 && i2 < *isize2)
 364     {
 365         /* prologue */
 366         rnr1 = atoms1->atom[index1[i1]].resind;
 367         rnr2 = atoms2->atom[index2[i2]].resind;
 368         if (rnr1 != prnr1 || rnr2 != prnr2)
 369         {
 370             if (debug)
 371             {
 372                 fprintf(debug, "R: %s%d %s%d\n",
 373                         *resinfo1[rnr1].name, rnr1, *resinfo2[rnr2].name, rnr2);
 374             }
 375             rescmp = std::strcmp(*resinfo1[rnr1].name, *resinfo2[rnr2].name);
 376         }
 377         if (debug)
 378         {
 379             fprintf(debug, "comparing %d %d", i1, i2);
 380         }
 381         atcmp = debug_strcmp(*atnms1[index1[i1]], *atnms2[index2[i2]]);
 382
 383         /* the works */
 384         if (atcmp != 0) /* no match -> find match within residues */
 385         {
 386             m1 = find_res_end(i1, *isize1, index1, atoms1);
 387             m2 = find_res_end(i2, *isize2, index2, atoms2);
 388             if (debug)
 389             {
 390                 fprintf(debug, " [%d<%d %d<%d]", i1, m1, i2, m2);
 391             }
 392             atcmp = find_next_match_atoms_in_res(&i1, index1, m1, atnms1,
 393                                                  &i2, index2, m2, atnms2);
 394             if (debug)
 395             {
 396                 fprintf(debug, " -> %d %d %s-%s", i1, i2,
 397                         *atnms1[index1[i1]], *atnms2[index2[i2]]);
 398             }
 399
 400         }
 401         if (atcmp != 0) /* still no match -> next residue(s) */
 402         {
 403             prnr1  = rnr1;
 404             prnr2  = rnr2;
 405             rescmp = find_next_match_res(&rnr1, rsize1, rindex1, resinfo1,
 406                                          &rnr2, rsize2, rindex2, resinfo2);
 407             if (rnr1 != prnr1)
 408             {
 409                 i1 = find_first_atom_in_res(rnr1, *isize1, index1, atoms1->atom);
 410             }
 411             if (rnr2 != prnr2)
 412             {
 413                 i2 = find_first_atom_in_res(rnr2, *isize2, index2, atoms2->atom);
 414             }
 415             if (debug)
 416             {
 417                 fprintf(debug, " -> %s%d-%s%d %s%d-%s%d",
 418                         *resinfo1[rnr1].name, rnr1,
 419                         *atnms1[index1[i1]], index1[i1],
 420                         *resinfo2[rnr2].name, rnr2,
 421                         *atnms2[index2[i2]], index2[i2]);
 422             }
 423             m1 = find_res_end(i1, *isize1, index1, atoms1);
 424             m2 = find_res_end(i2, *isize2, index2, atoms2);
 425             if (debug)
 426             {
 427                 fprintf(debug, " [%d<%d %d<%d]", i1, m1, i2, m2);
 428             }
 429             atcmp = find_next_match_atoms_in_res(&i1, index1, m1, atnms1,
 430                                                  &i2, index2, m2, atnms2);
 431             if (debug)
 432             {
 433                 fprintf(debug, " -> %d %d %s-%s", i1, i2,
 434                         *atnms1[index1[i1]], *atnms2[index2[i2]]);
 435             }
 436         }
 437         if (debug)
 438         {
 439             fprintf(debug, "(%d %d): %d %d\n", ii1, ii2, atcmp, rescmp);
 440         }
 441         if (atcmp == 0) /* if match -> store indices */
 442         {
 443             index1[ii1++] = index1[i1];
 444             index2[ii2++] = index2[i2];
 445         }
 446         i1++;
 447         i2++;
 448
 449         /* epilogue */
 450         prnr1 = rnr1;
 451         prnr2 = rnr2;
 452     }
 453
 454     if (ii1 != ii2)
 455     {
 456         gmx_fatal(FARGS, "DEATH HORROR: non-equal number of matching atoms!\n");
 457     }
 458     if (ii1 == i1 && ii2 == i2)
 459     {
 460         printf("All atoms in index groups 1 and 2 match\n");
 461     }
 462     else
 463     {
 464         if (i1 == i2 && ii1 == ii2)
 465         {
 466             printf("Both index groups modified from %d to %d atoms\n", i1, ii1);
 467         }
 468         else
 469         {
 470             if (ii1 != i1)
 471             {
 472                 printf("Index group 1 modified from %d to %d atoms\n", i1, ii1);
 473             }
 474             if (ii2 != i2)
 475             {
 476                 printf("Index group 2 modified from %d to %d atoms\n", i2, ii2);
 477             }
 478         }
 479         *isize1 = ii1;
 480         *isize2 = ii2;
 481     }
 482 }
 483 /* 1 */
 484
 485 int gmx_confrms(int argc, char *argv[])
 486 {
 487     const char     *desc[] = {
 488         "[THISMODULE] computes the root mean square deviation (RMSD) of two",
 489         "structures after least-squares fitting the second structure on the first one.",
 490         "The two structures do NOT need to have the same number of atoms,",
 491         "only the two index groups used for the fit need to be identical.",
 492         "With [TT]-name[tt] only matching atom names from the selected groups",
 493         "will be used for the fit and RMSD calculation. This can be useful ",
 494         "when comparing mutants of a protein.",
 495         "[PAR]",
 496         "The superimposed structures are written to file. In a [REF].pdb[ref] file",
 497         "the two structures will be written as separate models",
 498         "(use [TT]rasmol -nmrpdb[tt]). Also in a [REF].pdb[ref] file, B-factors",
 499         "calculated from the atomic MSD values can be written with [TT]-bfac[tt].",
 500     };
 501     static gmx_bool bOne  = FALSE, bRmpbc = FALSE, bMW = TRUE, bName = FALSE,
 502                     bBfac = FALSE, bFit = TRUE, bLabel = FALSE;
 503
 504     t_pargs  pa[] = {
 505         { "-one", FALSE, etBOOL, {&bOne},   "Only write the fitted structure to file" },
 506         { "-mw",  FALSE, etBOOL, {&bMW},    "Mass-weighted fitting and RMSD" },
 507         { "-pbc", FALSE, etBOOL, {&bRmpbc}, "Try to make molecules whole again" },
 508         { "-fit", FALSE, etBOOL, {&bFit},
 509           "Do least squares superposition of the target structure to the reference" },
 510         { "-name", FALSE, etBOOL, {&bName},
 511           "Only compare matching atom names" },
 512         { "-label", FALSE, etBOOL, {&bLabel},
 513           "Added chain labels A for first and B for second structure"},
 514         { "-bfac", FALSE, etBOOL, {&bBfac},
 515           "Output B-factors from atomic MSD values" }
 516     };
 517     t_filenm fnm[] = {
 518         { efTPS, "-f1",  "conf1.gro", ffREAD  },
 519         { efSTX, "-f2",  "conf2",     ffREAD  },
 520         { efSTO, "-o",   "fit.pdb",   ffWRITE },
 521         { efNDX, "-n1",  "fit1",      ffOPTRD },
 522         { efNDX, "-n2",  "fit2",      ffOPTRD },
 523         { efNDX, "-no",  "match",     ffOPTWR }
 524     };
 525 #define NFILE asize(fnm)
 526
 527     /* the two structure files */
 528     const char       *conf1file, *conf2file, *matchndxfile, *outfile;
 529     FILE             *fp;
 530     char             *name1, *name2;
 531     t_topology       *top1, *top2;
 532     int               ePBC1, ePBC2;
 533     t_atoms          *atoms1, *atoms2;
 534     int               warn = 0;
 535     int               at;
 536     real             *w_rls, mass, totmass;
 537     rvec             *x1, *v1, *x2, *v2, *fit_x;
 538     matrix            box1, box2;
 539
 540     gmx_output_env_t *oenv;
 541
 542     /* counters */
 543     int     i, m;
 544
 545     /* center of mass calculation */
 546     rvec    xcm1, xcm2;
 547
 548     /* variables for fit */
 549     char    *groupnames1, *groupnames2;
 550     int      isize1, isize2;
 551     int     *index1, *index2;
 552     real     rms, msd, minmsd, maxmsd;
 553     real    *msds;
 554
 555
 556     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW,
 557                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &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, top1, &ePBC1, &x1, &v1, box1, TRUE);
 569     atoms1 = &(top1->atoms);
 570     fprintf(stderr, "%s\nContaining %d atoms in %d residues\n",
 571             *top1->name, 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, top2, &ePBC2, &x2, &v2, box2, TRUE);
 592     atoms2 = &(top2->atoms);
 593     fprintf(stderr, "%s\nContaining %d atoms in %d residues\n",
 594             *top2->name, 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 = gmx_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, "%4d%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, "%4d%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 (std::strcmp(name1, name2))
 637         {
 638             if (warn < 20)
 639             {
 640                 fprintf(stderr,
 641                         "Warning: atomnames at index %d don't match: %d %s, %d %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 = nullptr;
 677     }
 678     else
 679     {
 680         clear_rvec(xcm1);
 681         clear_rvec(xcm2);
 682         w_rls = nullptr;
 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       = gmx::square(x1[index1[at]][m] - x2[index2[at]][m]);
 697             rms      += msd*mass;
 698             msds[at] += msd;
 699         }
 700         maxmsd   = std::max(maxmsd, msds[at]);
 701         minmsd   = std::min(minmsd, msds[at]);
 702         totmass += mass;
 703     }
 704     rms = std::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                 atoms1->havePdbInfo = TRUE;
 743
 744                 /* Avoid segfaults when writing the pdb-file */
 745                 for (i = 0; i < atoms1->nr; i++)
 746                 {
 747                     atoms1->pdbinfo[i].type         = eptAtom;
 748                     atoms1->pdbinfo[i].occup        = 1.00;
 749                     atoms1->pdbinfo[i].bAnisotropic = FALSE;
 750                     if (bBfac)
 751                     {
 752                         atoms1->pdbinfo[i].bfac = 0;
 753                     }
 754                     if (bLabel)
 755                     {
 756                         atoms1->resinfo[atoms1->atom[i].resind].chainid = 'A';
 757                     }
 758                 }
 759
 760                 for (i = 0; i < isize1; i++)
 761                 {
 762                     /* atoms1->pdbinfo[index1[i]].type = eptAtom; */
 763 /*  atoms1->pdbinfo[index1[i]].bAnisotropic = FALSE; */
 764                     if (bBfac)
 765                     {
 766                         atoms1->pdbinfo[index1[i]].bfac = (800*M_PI*M_PI/3.0)*msds[i];
 767                     }
 768 /*  if (bLabel) */
 769 /*    atoms1->resinfo[atoms1->atom[index1[i]].resind].chain = 'A'; */
 770                 }
 771                 srenew(atoms2->pdbinfo, atoms2->nr);
 772                 srenew(atoms2->atom, atoms2->nr); /* Why renew atom? */
 773
 774                 for (i = 0; i < atoms2->nr; i++)
 775                 {
 776                     atoms2->pdbinfo[i].type         = eptAtom;
 777                     atoms2->pdbinfo[i].occup        = 1.00;
 778                     atoms2->pdbinfo[i].bAnisotropic = FALSE;
 779                     if (bBfac)
 780                     {
 781                         atoms2->pdbinfo[i].bfac = 0;
 782                     }
 783                     if (bLabel)
 784                     {
 785                         atoms2->resinfo[atoms1->atom[i].resind].chainid = 'B';
 786                     }
 787                 }
 788
 789                 for (i = 0; i < isize2; i++)
 790                 {
 791                     /* atoms2->pdbinfo[index2[i]].type = eptAtom; */
 792 /*  atoms2->pdbinfo[index2[i]].bAnisotropic = FALSE; */
 793                     if (bBfac)
 794                     {
 795                         atoms2->pdbinfo[index2[i]].bfac = (800*M_PI*M_PI/3.0)*msds[i];
 796                     }
 797 /*  if (bLabel) */
 798 /*    atoms2->resinfo[atoms2->atom[index2[i]].resind].chain = 'B'; */
 799                 }
 800             }
 801             fp = gmx_ffopen(outfile, "w");
 802             if (!bOne)
 803             {
 804                 write_pdbfile(fp, *top1->name, atoms1, x1, ePBC1, box1, ' ', 1, nullptr, TRUE);
 805             }
 806             write_pdbfile(fp, *top2->name, atoms2, x2, ePBC2, box2, ' ', bOne ? -1 : 2, nullptr, TRUE);
 807             gmx_ffclose(fp);
 808             break;
 809         case efGRO:
 810             if (bBfac)
 811             {
 812                 fprintf(stderr, "WARNING: cannot write B-factor values to gro file\n");
 813             }
 814             fp = gmx_ffopen(outfile, "w");
 815             if (!bOne)
 816             {
 817                 write_hconf_p(fp, *top1->name, atoms1, x1, v1, box1);
 818             }
 819             write_hconf_p(fp, *top2->name, atoms2, x2, v2, box2);
 820             gmx_ffclose(fp);
 821             break;
 822         default:
 823             if (bBfac)
 824             {
 825                 fprintf(stderr, "WARNING: cannot write B-factor values to %s file\n",
 826                         ftp2ext(fn2ftp(outfile)));
 827             }
 828             if (!bOne)
 829             {
 830                 fprintf(stderr,
 831                         "WARNING: cannot write the reference structure to %s file\n",
 832                         ftp2ext(fn2ftp(outfile)));
 833             }
 834             write_sto_conf(outfile, *top2->name, atoms2, x2, v2, ePBC2, box2);
 835             break;
 836     }
 837
 838     view_all(oenv, NFILE, fnm);
 839
 840     return 0;
 841 }