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