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