src/tools/gmx_rms.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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * 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 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "smalloc.h"
  43 #include <math.h>
  44 #include "macros.h"
  45 #include "typedefs.h"
  46 #include "xvgr.h"
  47 #include "copyrite.h"
  48 #include "statutil.h"
  49 #include "string2.h"
  50 #include "vec.h"
  51 #include "index.h"
  52 #include "gmx_fatal.h"
  53 #include "futil.h"
  54 #include "princ.h"
  55 #include "rmpbc.h"
  56 #include "do_fit.h"
  57 #include "matio.h"
  58 #include "tpxio.h"
  59 #include "cmat.h"
  60 #include "viewit.h"
  61 #include "gmx_ana.h"
  62
  63
  64 static void norm_princ(t_atoms *atoms, int isize, atom_id *index, int natoms,
  65                        rvec *x)
  66 {
  67     int i, m;
  68     rvec princ, vec;
  69
  70     /* equalize principal components: */
  71     /* orient principal axes, get principal components */
  72     orient_princ(atoms, isize, index, natoms, x, NULL, princ);
  73
  74     /* calc our own principal components */
  75     clear_rvec(vec);
  76     for (m = 0; m < DIM; m++)
  77     {
  78         for (i = 0; i < isize; i++)
  79             vec[m] += sqr(x[index[i]][m]);
  80         vec[m] = sqrt(vec[m] / isize);
  81         /* calculate scaling constants */
  82         vec[m] = 1 / (sqrt(3) * vec[m]);
  83     }
  84
  85     /* scale coordinates */
  86     for (i = 0; i < natoms; i++)
  87     {
  88         for (m = 0; m < DIM; m++)
  89         {
  90             x[i][m] *= vec[m];
  91         }
  92     }
  93 }
  94
  95 int gmx_rms(int argc, char *argv[])
  96 {
  97     const char
  98         *desc[] =
  99             {
 100                 "[TT]g_rms[tt] compares two structures by computing the root mean square",
 101                 "deviation (RMSD), the size-independent [GRK]rho[grk] similarity parameter",
 102                 "([TT]rho[tt]) or the scaled [GRK]rho[grk] ([TT]rhosc[tt]), ",
 103                 "see Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).",
 104                 "This is selected by [TT]-what[tt].[PAR]"
 105
 106                 "Each structure from a trajectory ([TT]-f[tt]) is compared to a",
 107                 "reference structure. The reference structure",
 108                 "is taken from the structure file ([TT]-s[tt]).[PAR]",
 109
 110                 "With option [TT]-mir[tt] also a comparison with the mirror image of",
 111                 "the reference structure is calculated.",
 112                 "This is useful as a reference for 'significant' values, see",
 113                 "Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).[PAR]",
 114
 115                 "Option [TT]-prev[tt] produces the comparison with a previous frame",
 116                 "the specified number of frames ago.[PAR]",
 117
 118                 "Option [TT]-m[tt] produces a matrix in [TT].xpm[tt] format of",
 119                 "comparison values of each structure in the trajectory with respect to",
 120                 "each other structure. This file can be visualized with for instance",
 121                 "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].[PAR]",
 122
 123                 "Option [TT]-fit[tt] controls the least-squares fitting of",
 124                 "the structures on top of each other: complete fit (rotation and",
 125                 "translation), translation only, or no fitting at all.[PAR]",
 126
 127                 "Option [TT]-mw[tt] controls whether mass weighting is done or not.",
 128                 "If you select the option (default) and ",
 129                 "supply a valid [TT].tpr[tt] file masses will be taken from there, ",
 130                 "otherwise the masses will be deduced from the [TT]atommass.dat[tt] file in",
 131                 "[TT]GMXLIB[tt]. This is fine for proteins, but not",
 132                 "necessarily for other molecules. A default mass of 12.011 amu (carbon)",
 133                 "is assigned to unknown atoms. You can check whether this happend by",
 134                 "turning on the [TT]-debug[tt] flag and inspecting the log file.[PAR]",
 135
 136                 "With [TT]-f2[tt], the 'other structures' are taken from a second",
 137                 "trajectory, this generates a comparison matrix of one trajectory",
 138                 "versus the other.[PAR]",
 139
 140                 "Option [TT]-bin[tt] does a binary dump of the comparison matrix.[PAR]",
 141
 142                 "Option [TT]-bm[tt] produces a matrix of average bond angle deviations",
 143                 "analogously to the [TT]-m[tt] option. Only bonds between atoms in the",
 144                 "comparison group are considered." };
 145     static gmx_bool bPBC = TRUE, bFitAll = TRUE, bSplit = FALSE;
 146     static gmx_bool bDeltaLog = FALSE;
 147     static int prev = 0, freq = 1, freq2 = 1, nlevels = 80, avl = 0;
 148     static real rmsd_user_max = -1, rmsd_user_min = -1, bond_user_max = -1,
 149         bond_user_min = -1, delta_maxy = 0.0;
 150     /* strings and things for selecting difference method */
 151     enum
 152     {
 153         ewSel, ewRMSD, ewRho, ewRhoSc, ewNR
 154     };
 155     int ewhat;
 156     const char *what[ewNR + 1] =
 157         { NULL, "rmsd", "rho", "rhosc", NULL };
 158     const char *whatname[ewNR] =
 159         { NULL, "RMSD", "Rho", "Rho sc" };
 160     const char *whatlabel[ewNR] =
 161         { NULL, "RMSD (nm)", "Rho", "Rho sc" };
 162     const char *whatxvgname[ewNR] =
 163         { NULL, "RMSD", "\\8r\\4", "\\8r\\4\\ssc\\N" };
 164     const char *whatxvglabel[ewNR] =
 165         { NULL, "RMSD (nm)", "\\8r\\4", "\\8r\\4\\ssc\\N" };
 166     /* strings and things for fitting methods */
 167     enum
 168     {
 169         efSel, efFit, efReset, efNone, efNR
 170     };
 171     int efit;
 172     const char *fit[efNR + 1] =
 173         { NULL, "rot+trans", "translation", "none", NULL };
 174     const char *fitgraphlabel[efNR + 1] =
 175         { NULL, "lsq fit", "translational fit", "no fit" };
 176     static int nrms = 1;
 177     static gmx_bool bMassWeighted = TRUE;
 178     t_pargs pa[] =
 179         {
 180             { "-what", FALSE, etENUM,
 181                 { what }, "Structural difference measure" },
 182             { "-pbc", FALSE, etBOOL,
 183                 { &bPBC }, "PBC check" },
 184             { "-fit", FALSE, etENUM,
 185                 { fit }, "Fit to reference structure" },
 186             { "-prev", FALSE, etINT,
 187                 { &prev }, "Compare with previous frame" },
 188             { "-split", FALSE, etBOOL,
 189                 { &bSplit }, "Split graph where time is zero" },
 190             { "-fitall", FALSE, etBOOL,
 191                 { &bFitAll }, "HIDDENFit all pairs of structures in matrix" },
 192             { "-skip", FALSE, etINT,
 193                 { &freq }, "Only write every nr-th frame to matrix" },
 194             { "-skip2", FALSE, etINT,
 195                 { &freq2 }, "Only write every nr-th frame to matrix" },
 196             { "-max", FALSE, etREAL,
 197                 { &rmsd_user_max }, "Maximum level in comparison matrix" },
 198             { "-min", FALSE, etREAL,
 199                 { &rmsd_user_min }, "Minimum level in comparison matrix" },
 200             { "-bmax", FALSE, etREAL,
 201                 { &bond_user_max }, "Maximum level in bond angle matrix" },
 202             { "-bmin", FALSE, etREAL,
 203                 { &bond_user_min }, "Minimum level in bond angle matrix" },
 204             { "-mw", FALSE, etBOOL,
 205                 { &bMassWeighted }, "Use mass weighting for superposition" },
 206             { "-nlevels", FALSE, etINT,
 207                 { &nlevels }, "Number of levels in the matrices" },
 208             { "-ng", FALSE, etINT,
 209                 { &nrms }, "Number of groups to compute RMS between" },
 210             { "-dlog", FALSE, etBOOL,
 211                 { &bDeltaLog },
 212                 "HIDDENUse a log x-axis in the delta t matrix" },
 213             { "-dmax", FALSE, etREAL,
 214                 { &delta_maxy }, "HIDDENMaximum level in delta matrix" },
 215             { "-aver", FALSE, etINT,
 216                 { &avl },
 217                 "HIDDENAverage over this distance in the RMSD matrix" } };
 218     int natoms_trx, natoms_trx2, natoms;
 219     int i, j, k, m, teller, teller2, tel_mat, tel_mat2;
 220 #define NFRAME 5000
 221     int maxframe = NFRAME, maxframe2 = NFRAME;
 222     real t, *w_rls, *w_rms, *w_rls_m = NULL, *w_rms_m = NULL;
 223     gmx_bool bNorm, bAv, bFreq2, bFile2, bMat, bBond, bDelta, bMirror, bMass;
 224     gmx_bool bFit, bReset;
 225     t_topology top;
 226     int ePBC;
 227     t_iatom *iatom = NULL;
 228
 229     matrix box;
 230     rvec *x, *xp, *xm = NULL, **mat_x = NULL, **mat_x2, *mat_x2_j = NULL, vec1,
 231         vec2;
 232     t_trxstatus *status;
 233     char buf[256], buf2[256];
 234     int ncons = 0;
 235     FILE *fp;
 236     real rlstot = 0, **rls, **rlsm = NULL, *time, *time2, *rlsnorm = NULL,
 237         **rmsd_mat = NULL, **bond_mat = NULL, *axis, *axis2, *del_xaxis,
 238         *del_yaxis, rmsd_max, rmsd_min, rmsd_avg, bond_max, bond_min, ang;
 239     real **rmsdav_mat = NULL, av_tot, weight, weight_tot;
 240     real **delta = NULL, delta_max, delta_scalex = 0, delta_scaley = 0,
 241         *delta_tot;
 242     int delta_xsize = 0, del_lev = 100, mx, my, abs_my;
 243     gmx_bool bA1, bA2, bPrev, bTop, *bInMat = NULL;
 244     int ifit, *irms, ibond = 0, *ind_bond1 = NULL, *ind_bond2 = NULL, n_ind_m =
 245         0;
 246     atom_id *ind_fit, **ind_rms, *ind_m = NULL, *rev_ind_m = NULL, *ind_rms_m =
 247         NULL;
 248     char *gn_fit, **gn_rms;
 249     t_rgb rlo, rhi;
 250     output_env_t oenv;
 251     gmx_rmpbc_t  gpbc=NULL;
 252
 253     t_filenm fnm[] =
 254         {
 255             { efTPS, NULL, NULL, ffREAD },
 256             { efTRX, "-f", NULL, ffREAD },
 257             { efTRX, "-f2", NULL, ffOPTRD },
 258             { efNDX, NULL, NULL, ffOPTRD },
 259             { efXVG, NULL, "rmsd", ffWRITE },
 260             { efXVG, "-mir", "rmsdmir", ffOPTWR },
 261             { efXVG, "-a", "avgrp", ffOPTWR },
 262             { efXVG, "-dist", "rmsd-dist", ffOPTWR },
 263             { efXPM, "-m", "rmsd", ffOPTWR },
 264             { efDAT, "-bin", "rmsd", ffOPTWR },
 265             { efXPM, "-bm", "bond", ffOPTWR } };
 266 #define NFILE asize(fnm)
 267
 268     CopyRight(stderr, argv[0]);
 269     parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW
 270         | PCA_BE_NICE, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
 271                       &oenv);
 272     /* parse enumerated options: */
 273     ewhat = nenum(what);
 274     if (ewhat == ewRho || ewhat == ewRhoSc)
 275         please_cite(stdout, "Maiorov95");
 276     efit = nenum(fit);
 277     bFit = efit == efFit;
 278     bReset = efit == efReset;
 279     if (bFit)
 280     {
 281         bReset = TRUE; /* for fit, reset *must* be set */
 282     }
 283     else
 284     {
 285         bFitAll = FALSE;
 286     }
 287
 288     /* mark active cmdline options */
 289     bMirror = opt2bSet("-mir", NFILE, fnm); /* calc RMSD vs mirror of ref. */
 290     bFile2 = opt2bSet("-f2", NFILE, fnm);
 291     bMat = opt2bSet("-m", NFILE, fnm);
 292     bBond = opt2bSet("-bm", NFILE, fnm);
 293     bDelta = (delta_maxy > 0); /* calculate rmsd vs delta t matrix from *
 294      *  your RMSD matrix (hidden option       */
 295     bNorm = opt2bSet("-a", NFILE, fnm);
 296     bFreq2 = opt2parg_bSet("-skip2", asize(pa), pa);
 297     if (freq <= 0)
 298     {
 299         fprintf(stderr, "The number of frames to skip is <= 0. "
 300             "Writing out all frames.\n\n");
 301         freq = 1;
 302     }
 303     if (!bFreq2)
 304     {
 305         freq2 = freq;
 306     }
 307     else if (bFile2 && freq2 <= 0)
 308     {
 309         fprintf(stderr,
 310                 "The number of frames to skip in second trajectory is <= 0.\n"
 311                     "  Writing out all frames.\n\n");
 312         freq2 = 1;
 313     }
 314
 315     bPrev = (prev > 0);
 316     if (bPrev)
 317     {
 318         prev = abs(prev);
 319         if (freq != 1)
 320             fprintf(stderr, "WARNING: option -skip also applies to -prev\n");
 321     }
 322
 323     if (bFile2 && !bMat && !bBond)
 324     {
 325         fprintf(
 326                 stderr,
 327                 "WARNING: second trajectory (-f2) useless when not calculating matrix (-m/-bm),\n"
 328                     "         will not read from %s\n", opt2fn("-f2", NFILE,
 329                                                                fnm));
 330         bFile2 = FALSE;
 331     }
 332
 333     if (bDelta)
 334     {
 335         bMat = TRUE;
 336         if (bFile2)
 337         {
 338             fprintf(stderr,
 339                     "WARNING: second trajectory (-f2) useless when making delta matrix,\n"
 340                         "         will not read from %s\n", opt2fn("-f2",
 341                                                                    NFILE, fnm));
 342             bFile2 = FALSE;
 343         }
 344     }
 345
 346     bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xp,
 347                          NULL, box, TRUE);
 348     snew(w_rls,top.atoms.nr);
 349     snew(w_rms,top.atoms.nr);
 350
 351     if (!bTop && bBond)
 352     {
 353         fprintf(stderr,
 354                 "WARNING: Need a run input file for bond angle matrix,\n"
 355                     "         will not calculate bond angle matrix.\n");
 356         bBond = FALSE;
 357     }
 358
 359     if (bReset)
 360     {
 361         fprintf(stderr, "Select group for %s fit\n", bFit ? "least squares"
 362             : "translational");
 363         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &ifit,
 364                   &ind_fit, &gn_fit);
 365     }
 366     else
 367         ifit = 0;
 368
 369     if (bReset)
 370     {
 371         if (bFit && ifit < 3)
 372             gmx_fatal(FARGS,"Need >= 3 points to fit!\n" );
 373
 374         bMass = FALSE;
 375         for(i=0; i<ifit; i++)
 376         {
 377             if (bMassWeighted)
 378                 w_rls[ind_fit[i]] = top.atoms.atom[ind_fit[i]].m;
 379             else
 380                 w_rls[ind_fit[i]] = 1;
 381             bMass = bMass || (top.atoms.atom[ind_fit[i]].m != 0);
 382         }
 383         if (!bMass)
 384         {
 385             fprintf(stderr,"All masses in the fit group are 0, using masses of 1\n");
 386             for(i=0; i<ifit; i++)
 387             {
 388                 w_rls[ind_fit[i]] = 1;
 389             }
 390         }
 391     }
 392
 393     if (bMat || bBond)
 394         nrms=1;
 395
 396     snew(gn_rms,nrms);
 397     snew(ind_rms,nrms);
 398     snew(irms,nrms);
 399
 400     fprintf(stderr,"Select group%s for %s calculation\n",
 401             (nrms>1) ? "s" : "",whatname[ewhat]);
 402     get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
 403               nrms,irms,ind_rms,gn_rms);
 404
 405     if (bNorm)
 406     {
 407         snew(rlsnorm,irms[0]);
 408     }
 409     snew(rls,nrms);
 410     for(j=0; j<nrms; j++)
 411         snew(rls[j],maxframe);
 412     if (bMirror)
 413     {
 414         snew(rlsm,nrms);
 415         for(j=0; j<nrms; j++)
 416             snew(rlsm[j],maxframe);
 417     }
 418     snew(time,maxframe);
 419     for(j=0; j<nrms; j++)
 420     {
 421         bMass = FALSE;
 422         for(i=0; i<irms[j]; i++)
 423         {
 424             if (bMassWeighted)
 425                 w_rms[ind_rms[j][i]] = top.atoms.atom[ind_rms[j][i]].m;
 426             else
 427                 w_rms[ind_rms[j][i]] = 1.0;
 428             bMass = bMass || (top.atoms.atom[ind_rms[j][i]].m != 0);
 429         }
 430         if (!bMass) {
 431             fprintf(stderr,"All masses in group %d are 0, using masses of 1\n",j);
 432             for(i=0; i<irms[j]; i++)
 433                 w_rms[ind_rms[j][i]] = 1;
 434         }
 435     }
 436     /* Prepare reference frame */
 437     if (bPBC) {
 438       gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,box);
 439       gmx_rmpbc(gpbc,top.atoms.nr,box,xp);
 440     }
 441     if (bReset)
 442         reset_x(ifit,ind_fit,top.atoms.nr,NULL,xp,w_rls);
 443     if (bMirror) {
 444         /* generate reference structure mirror image: */
 445         snew(xm, top.atoms.nr);
 446         for(i=0; i<top.atoms.nr; i++) {
 447             copy_rvec(xp[i],xm[i]);
 448             xm[i][XX] = -xm[i][XX];
 449         }
 450     }
 451     if (ewhat==ewRhoSc)
 452         norm_princ(&top.atoms, ifit, ind_fit, top.atoms.nr, xp);
 453
 454     /* read first frame */
 455     natoms_trx=read_first_x(oenv,&status,opt2fn("-f",NFILE,fnm),&t,&x,box);
 456     if (natoms_trx != top.atoms.nr)
 457         fprintf(stderr,
 458                 "\nWARNING: topology has %d atoms, whereas trajectory has %d\n",
 459                 top.atoms.nr,natoms_trx);
 460     natoms = min(top.atoms.nr,natoms_trx);
 461     if (bMat || bBond || bPrev) {
 462         snew(mat_x,NFRAME);
 463
 464         if (bPrev)
 465             /* With -prev we use all atoms for simplicity */
 466             n_ind_m = natoms;
 467         else {
 468             /* Check which atoms we need (fit/rms) */
 469             snew(bInMat,natoms);
 470             for(i=0; i<ifit; i++)
 471                 bInMat[ind_fit[i]] = TRUE;
 472             n_ind_m=ifit;
 473             for(i=0; i<irms[0]; i++)
 474                 if (!bInMat[ind_rms[0][i]]) {
 475                     bInMat[ind_rms[0][i]] = TRUE;
 476                     n_ind_m++;
 477                 }
 478         }
 479         /* Make an index of needed atoms */
 480         snew(ind_m,n_ind_m);
 481         snew(rev_ind_m,natoms);
 482         j = 0;
 483         for(i=0; i<natoms; i++)
 484             if (bPrev || bInMat[i]) {
 485                 ind_m[j] = i;
 486                 rev_ind_m[i] = j;
 487                 j++;
 488             }
 489         snew(w_rls_m,n_ind_m);
 490         snew(ind_rms_m,irms[0]);
 491         snew(w_rms_m,n_ind_m);
 492         for(i=0; i<ifit; i++)
 493             w_rls_m[rev_ind_m[ind_fit[i]]] = w_rls[ind_fit[i]];
 494         for(i=0; i<irms[0]; i++) {
 495             ind_rms_m[i] = rev_ind_m[ind_rms[0][i]];
 496             w_rms_m[ind_rms_m[i]] = w_rms[ind_rms[0][i]];
 497         }
 498         sfree(bInMat);
 499     }
 500
 501     if (bBond) {
 502         ncons = 0;
 503         for(k=0; k<F_NRE; k++)
 504             if (IS_CHEMBOND(k)) {
 505                 iatom  = top.idef.il[k].iatoms;
 506                 ncons += top.idef.il[k].nr/3;
 507             }
 508         fprintf(stderr,"Found %d bonds in topology\n",ncons);
 509         snew(ind_bond1,ncons);
 510         snew(ind_bond2,ncons);
 511         ibond=0;
 512         for(k=0; k<F_NRE; k++)
 513             if (IS_CHEMBOND(k)) {
 514                 iatom = top.idef.il[k].iatoms;
 515                 ncons = top.idef.il[k].nr/3;
 516                 for (i=0; i<ncons; i++) {
 517                     bA1=FALSE;
 518                     bA2=FALSE;
 519                     for (j=0; j<irms[0]; j++) {
 520                         if (iatom[3*i+1] == ind_rms[0][j]) bA1=TRUE;
 521                         if (iatom[3*i+2] == ind_rms[0][j]) bA2=TRUE;
 522                     }
 523                     if (bA1 && bA2) {
 524                         ind_bond1[ibond] = rev_ind_m[iatom[3*i+1]];
 525                         ind_bond2[ibond] = rev_ind_m[iatom[3*i+2]];
 526                         ibond++;
 527                     }
 528                 }
 529             }
 530         fprintf(stderr,"Using %d bonds for bond angle matrix\n",ibond);
 531         if (ibond==0)
 532             gmx_fatal(FARGS,"0 bonds found");
 533     }
 534
 535     /* start looping over frames: */
 536     tel_mat = 0;
 537     teller = 0;
 538     do {
 539         if (bPBC)
 540           gmx_rmpbc(gpbc,natoms,box,x);
 541
 542         if (bReset)
 543             reset_x(ifit,ind_fit,natoms,NULL,x,w_rls);
 544         if (ewhat==ewRhoSc)
 545             norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 546
 547         if (bFit)
 548             /*do the least squares fit to original structure*/
 549             do_fit(natoms,w_rls,xp,x);
 550
 551         if (teller % freq == 0) {
 552             /* keep frame for matrix calculation */
 553             if (bMat || bBond || bPrev) {
 554                 if (tel_mat >= NFRAME)
 555                     srenew(mat_x,tel_mat+1);
 556                 snew(mat_x[tel_mat],n_ind_m);
 557                 for (i=0;i<n_ind_m;i++)
 558                     copy_rvec(x[ind_m[i]],mat_x[tel_mat][i]);
 559             }
 560             tel_mat++;
 561         }
 562
 563         /*calculate energy of root_least_squares*/
 564         if (bPrev) {
 565             j=tel_mat-prev-1;
 566             if (j<0)
 567                 j=0;
 568             for (i=0;i<n_ind_m;i++)
 569                 copy_rvec(mat_x[j][i],xp[ind_m[i]]);
 570             if (bReset)
 571                 reset_x(ifit,ind_fit,natoms,NULL,xp,w_rls);
 572             if (bFit)
 573                 do_fit(natoms,w_rls,x,xp);
 574         }
 575         for(j=0; (j<nrms); j++)
 576             rls[j][teller] =
 577                 calc_similar_ind(ewhat!=ewRMSD,irms[j],ind_rms[j],w_rms,x,xp);
 578         if (bNorm) {
 579             for(j=0; (j<irms[0]); j++)
 580                 rlsnorm[j] +=
 581                     calc_similar_ind(ewhat!=ewRMSD,1,&(ind_rms[0][j]),w_rms,x,xp);
 582         }
 583
 584         if (bMirror) {
 585             if (bFit)
 586                 /*do the least squares fit to mirror of original structure*/
 587                 do_fit(natoms,w_rls,xm,x);
 588
 589             for(j=0; j<nrms; j++)
 590                 rlsm[j][teller] =
 591                     calc_similar_ind(ewhat!=ewRMSD,irms[j],ind_rms[j],w_rms,x,xm);
 592         }
 593         time[teller]=output_env_conv_time(oenv,t);
 594
 595         teller++;
 596         if (teller >= maxframe) {
 597             maxframe +=NFRAME;
 598             srenew(time,maxframe);
 599             for(j=0; (j<nrms); j++)
 600                 srenew(rls[j],maxframe);
 601             if (bMirror)
 602                 for(j=0; (j<nrms); j++)
 603                     srenew(rlsm[j],maxframe);
 604         }
 605     } while (read_next_x(oenv,status,&t,natoms_trx,x,box));
 606     close_trj(status);
 607
 608     if (bFile2) {
 609         snew(time2,maxframe2);
 610
 611         fprintf(stderr,"\nWill read second trajectory file\n");
 612         snew(mat_x2,NFRAME);
 613         natoms_trx2 =
 614           read_first_x(oenv,&status,opt2fn("-f2",NFILE,fnm),&t,&x,box);
 615         if ( natoms_trx2 != natoms_trx )
 616             gmx_fatal(FARGS,
 617                       "Second trajectory (%d atoms) does not match the first one"
 618                       " (%d atoms)", natoms_trx2, natoms_trx);
 619         tel_mat2 = 0;
 620         teller2 = 0;
 621         do {
 622             if (bPBC)
 623               gmx_rmpbc(gpbc,natoms,box,x);
 624
 625             if (bReset)
 626                 reset_x(ifit,ind_fit,natoms,NULL,x,w_rls);
 627             if (ewhat==ewRhoSc)
 628                 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 629
 630             if (bFit)
 631                 /*do the least squares fit to original structure*/
 632                 do_fit(natoms,w_rls,xp,x);
 633
 634             if (teller2 % freq2 == 0) {
 635                 /* keep frame for matrix calculation */
 636                 if (bMat) {
 637                     if (tel_mat2 >= NFRAME)
 638                         srenew(mat_x2,tel_mat2+1);
 639                     snew(mat_x2[tel_mat2],n_ind_m);
 640                     for (i=0;i<n_ind_m;i++)
 641                         copy_rvec(x[ind_m[i]],mat_x2[tel_mat2][i]);
 642                 }
 643                 tel_mat2++;
 644             }
 645
 646             time2[teller2]=output_env_conv_time(oenv,t);
 647
 648             teller2++;
 649             if (teller2 >= maxframe2) {
 650                 maxframe2 +=NFRAME;
 651                 srenew(time2,maxframe2);
 652             }
 653         } while (read_next_x(oenv,status,&t,natoms_trx2,x,box));
 654         close_trj(status);
 655     } else {
 656         mat_x2=mat_x;
 657         time2=time;
 658         tel_mat2=tel_mat;
 659         freq2=freq;
 660     }
 661     gmx_rmpbc_done(gpbc);
 662
 663     if (bMat || bBond) {
 664         /* calculate RMS matrix */
 665         fprintf(stderr,"\n");
 666         if (bMat) {
 667             fprintf(stderr,"Building %s matrix, %dx%d elements\n",
 668                     whatname[ewhat],tel_mat,tel_mat2);
 669             snew(rmsd_mat,tel_mat);
 670         }
 671         if (bBond) {
 672             fprintf(stderr,"Building bond angle matrix, %dx%d elements\n",
 673                     tel_mat,tel_mat2);
 674             snew(bond_mat,tel_mat);
 675         }
 676         snew(axis,tel_mat);
 677         snew(axis2,tel_mat2);
 678         rmsd_max=0;
 679         if (bFile2)
 680             rmsd_min=1e10;
 681         else
 682             rmsd_min=0;
 683         rmsd_avg=0;
 684         bond_max=0;
 685         bond_min=1e10;
 686         for(j=0; j<tel_mat2; j++)
 687             axis2[j]=time2[freq2*j];
 688         if (bDelta) {
 689             if (bDeltaLog) {
 690                 delta_scalex=8.0/log(2.0);
 691                 delta_xsize=(int)(log(tel_mat/2)*delta_scalex+0.5)+1;
 692             }
 693             else {
 694                 delta_xsize=tel_mat/2;
 695             }
 696             delta_scaley=1.0/delta_maxy;
 697             snew(delta,delta_xsize);
 698             for(j=0; j<delta_xsize; j++)
 699                 snew(delta[j],del_lev+1);
 700             if (avl > 0) {
 701                 snew(rmsdav_mat,tel_mat);
 702                 for(j=0; j<tel_mat; j++)
 703                     snew(rmsdav_mat[j],tel_mat);
 704             }
 705         }
 706
 707         if (bFitAll)
 708             snew(mat_x2_j,natoms);
 709         for(i=0; i<tel_mat; i++) {
 710             axis[i]=time[freq*i];
 711             fprintf(stderr,"\r element %5d; time %5.2f  ",i,axis[i]);
 712             if (bMat) snew(rmsd_mat[i],tel_mat2);
 713             if (bBond) snew(bond_mat[i],tel_mat2);
 714             for(j=0; j<tel_mat2; j++) {
 715                 if (bFitAll) {
 716                     for (k=0;k<n_ind_m;k++)
 717                         copy_rvec(mat_x2[j][k],mat_x2_j[k]);
 718                     do_fit(n_ind_m,w_rls_m,mat_x[i],mat_x2_j);
 719                 } else
 720                     mat_x2_j=mat_x2[j];
 721                 if (bMat) {
 722                     if (bFile2 || (i<j)) {
 723                         rmsd_mat[i][j] =
 724                             calc_similar_ind(ewhat!=ewRMSD,irms[0],ind_rms_m,
 725                                              w_rms_m,mat_x[i],mat_x2_j);
 726                         if (rmsd_mat[i][j] > rmsd_max) rmsd_max=rmsd_mat[i][j];
 727                         if (rmsd_mat[i][j] < rmsd_min) rmsd_min=rmsd_mat[i][j];
 728                         rmsd_avg += rmsd_mat[i][j];
 729                     } else
 730                         rmsd_mat[i][j]=rmsd_mat[j][i];
 731                 }
 732                 if (bBond) {
 733                     if (bFile2 || (i<=j)) {
 734                         ang=0.0;
 735                         for(m=0;m<ibond;m++) {
 736                             rvec_sub(mat_x[i][ind_bond1[m]],mat_x[i][ind_bond2[m]],vec1);
 737                             rvec_sub(mat_x2_j[ind_bond1[m]],mat_x2_j[ind_bond2[m]],vec2);
 738                             ang += acos(cos_angle(vec1,vec2));
 739                         }
 740                         bond_mat[i][j]=ang*180.0/(M_PI*ibond);
 741                         if (bond_mat[i][j] > bond_max) bond_max=bond_mat[i][j];
 742                         if (bond_mat[i][j] < bond_min) bond_min=bond_mat[i][j];
 743                     }
 744                     else
 745                         bond_mat[i][j]=bond_mat[j][i];
 746                 }
 747             }
 748         }
 749         if (bFile2)
 750             rmsd_avg /= tel_mat*tel_mat2;
 751         else
 752             rmsd_avg /= tel_mat*(tel_mat - 1)/2;
 753         if (bMat && (avl > 0)) {
 754             rmsd_max=0.0;
 755             rmsd_min=0.0;
 756             rmsd_avg=0.0;
 757             for(j=0; j<tel_mat-1; j++) {
 758                 for(i=j+1; i<tel_mat; i++) {
 759                     av_tot=0;
 760                     weight_tot=0;
 761                     for (my=-avl; my<=avl; my++) {
 762                         if ((j+my>=0) && (j+my<tel_mat)) {
 763                             abs_my = abs(my);
 764                             for (mx=-avl; mx<=avl; mx++) {
 765                                 if ((i+mx>=0) && (i+mx<tel_mat)) {
 766                                     weight = (real)(avl+1-max(abs(mx),abs_my));
 767                                     av_tot += weight*rmsd_mat[i+mx][j+my];
 768                                     weight_tot+=weight;
 769                                 }
 770                             }
 771                         }
 772                     }
 773                     rmsdav_mat[i][j] = av_tot/weight_tot;
 774                     rmsdav_mat[j][i] = rmsdav_mat[i][j];
 775                     if (rmsdav_mat[i][j] > rmsd_max) rmsd_max=rmsdav_mat[i][j];
 776                 }
 777             }
 778             rmsd_mat=rmsdav_mat;
 779         }
 780
 781         if (bMat) {
 782             fprintf(stderr,"\n%s: Min %f, Max %f, Avg %f\n",
 783                     whatname[ewhat],rmsd_min,rmsd_max,rmsd_avg);
 784             rlo.r = 1; rlo.g = 1; rlo.b = 1;
 785             rhi.r = 0; rhi.g = 0; rhi.b = 0;
 786             if (rmsd_user_max != -1) rmsd_max=rmsd_user_max;
 787             if (rmsd_user_min != -1) rmsd_min=rmsd_user_min;
 788             if ((rmsd_user_max !=  -1) || (rmsd_user_min != -1))
 789                 fprintf(stderr,"Min and Max value set to resp. %f and %f\n",
 790                         rmsd_min,rmsd_max);
 791             sprintf(buf,"%s %s matrix",gn_rms[0],whatname[ewhat]);
 792             write_xpm(opt2FILE("-m",NFILE,fnm,"w"),0,buf,whatlabel[ewhat],
 793                       output_env_get_time_label(oenv),output_env_get_time_label(oenv),tel_mat,tel_mat2,
 794                       axis,axis2, rmsd_mat,rmsd_min,rmsd_max,rlo,rhi,&nlevels);
 795             /* Print the distribution of RMSD values */
 796             if (opt2bSet("-dist",NFILE,fnm))
 797                 low_rmsd_dist(opt2fn("-dist",NFILE,fnm),rmsd_max,tel_mat,rmsd_mat,oenv);
 798
 799             if (bDelta) {
 800                 snew(delta_tot,delta_xsize);
 801                 for(j=0; j<tel_mat-1; j++) {
 802                     for(i=j+1; i<tel_mat; i++) {
 803                         mx=i-j ;
 804                         if (mx < tel_mat/2) {
 805                             if (bDeltaLog)
 806                                 mx=(int)(log(mx)*delta_scalex+0.5);
 807                             my=(int)(rmsd_mat[i][j]*delta_scaley*del_lev+0.5);
 808                             delta_tot[mx] += 1.0;
 809                             if ((rmsd_mat[i][j]>=0) && (rmsd_mat[i][j]<=delta_maxy))
 810                                 delta[mx][my] += 1.0;
 811                         }
 812                     }
 813                 }
 814                 delta_max=0;
 815                 for(i=0; i<delta_xsize; i++) {
 816                     if (delta_tot[i] > 0.0) {
 817                         delta_tot[i]=1.0/delta_tot[i];
 818                         for(j=0; j<=del_lev; j++) {
 819                             delta[i][j] *= delta_tot[i];
 820                             if (delta[i][j] > delta_max)
 821                                 delta_max=delta[i][j];
 822                         }
 823                     }
 824                 }
 825                 fprintf(stderr,"Maximum in delta matrix: %f\n",delta_max);
 826                 snew(del_xaxis,delta_xsize);
 827                 snew(del_yaxis,del_lev+1);
 828                 for (i=0; i<delta_xsize; i++)
 829                     del_xaxis[i]=axis[i]-axis[0];
 830                 for (i=0; i<del_lev+1; i++)
 831                     del_yaxis[i]=delta_maxy*i/del_lev;
 832                 sprintf(buf,"%s %s vs. delta t",gn_rms[0],whatname[ewhat]);
 833                 fp = ffopen("delta.xpm","w");
 834                 write_xpm(fp,0,buf,"density",output_env_get_time_label(oenv),whatlabel[ewhat],
 835                           delta_xsize,del_lev+1,del_xaxis,del_yaxis,
 836                           delta,0.0,delta_max,rlo,rhi,&nlevels);
 837                 ffclose(fp);
 838             }
 839             if (opt2bSet("-bin",NFILE,fnm)) {
 840                 /* NB: File must be binary if we use fwrite */
 841                 fp=ftp2FILE(efDAT,NFILE,fnm,"wb");
 842                 for(i=0;i<tel_mat;i++)
 843                     if(fwrite(rmsd_mat[i],sizeof(**rmsd_mat),tel_mat2,fp) != tel_mat2)
 844                     {
 845                         gmx_fatal(FARGS,"Error writing to output file");
 846                     }
 847                 ffclose(fp);
 848             }
 849         }
 850         if (bBond) {
 851             fprintf(stderr,"\nMin. angle: %f, Max. angle: %f\n",bond_min,bond_max);
 852             if (bond_user_max != -1) bond_max=bond_user_max;
 853             if (bond_user_min != -1) bond_min=bond_user_min;
 854             if ((bond_user_max !=  -1) || (bond_user_min != -1))
 855                 fprintf(stderr,"Bond angle Min and Max set to:\n"
 856                         "Min. angle: %f, Max. angle: %f\n",bond_min,bond_max);
 857             rlo.r = 1; rlo.g = 1; rlo.b = 1;
 858             rhi.r = 0; rhi.g = 0; rhi.b = 0;
 859             sprintf(buf,"%s av. bond angle deviation",gn_rms[0]);
 860             write_xpm(opt2FILE("-bm",NFILE,fnm,"w"),0,buf,"degrees",
 861                       output_env_get_time_label(oenv),output_env_get_time_label(oenv),tel_mat,tel_mat2,
 862                       axis,axis2, bond_mat,bond_min,bond_max,rlo,rhi,&nlevels);
 863         }
 864     }
 865
 866     bAv=opt2bSet("-a",NFILE,fnm);
 867
 868     /* Write the RMSD's to file */
 869     if (!bPrev)
 870         sprintf(buf,"%s",whatxvgname[ewhat]);
 871     else
 872         sprintf(buf,"%s with frame %g %s ago",whatxvgname[ewhat],
 873                 time[prev*freq]-time[0], output_env_get_time_label(oenv));
 874     fp=xvgropen(opt2fn("-o",NFILE,fnm),buf,output_env_get_xvgr_tlabel(oenv),
 875                 whatxvglabel[ewhat], oenv);
 876     if (output_env_get_print_xvgr_codes(oenv))
 877         fprintf(fp,"@ subtitle \"%s%s after %s%s%s\"\n",
 878                 (nrms==1)?"":"of "    , gn_rms[0], fitgraphlabel[efit],
 879                     bFit     ?" to ":""   , bFit?gn_fit:"");
 880     if (nrms != 1)
 881         xvgr_legend(fp,nrms,(const char**)gn_rms,oenv);
 882     for(i=0; (i<teller); i++) {
 883         if ( bSplit && i>0 &&
 884             abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv))<1e-5 )
 885         {
 886             fprintf(fp,"&\n");
 887         }
 888         fprintf(fp,"%12.7f",time[bPrev ? freq*i : i]);
 889         for(j=0; (j<nrms); j++) {
 890             fprintf(fp," %12.7f",rls[j][i]);
 891             if (bAv)
 892                 rlstot+=rls[j][i];
 893         }
 894         fprintf(fp,"\n");
 895     }
 896     ffclose(fp);
 897
 898     if (bMirror) {
 899         /* Write the mirror RMSD's to file */
 900         sprintf(buf,"%s with Mirror",whatxvgname[ewhat]);
 901         sprintf(buf2,"Mirror %s",whatxvglabel[ewhat]);
 902         fp=xvgropen(opt2fn("-mir",NFILE,fnm), buf, output_env_get_xvgr_tlabel(oenv),
 903                     buf2,oenv);
 904         if (nrms == 1) {
 905             if (output_env_get_print_xvgr_codes(oenv))
 906                 fprintf(fp,"@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
 907                         gn_rms[0],gn_fit);
 908         }
 909         else {
 910             if (output_env_get_print_xvgr_codes(oenv))
 911                 fprintf(fp,"@ subtitle \"after lsq fit to mirror %s\"\n",gn_fit);
 912             xvgr_legend(fp,nrms,(const char**)gn_rms,oenv);
 913         }
 914         for(i=0; (i<teller); i++) {
 915             if ( bSplit && i>0 && abs(time[i])<1e-5 )
 916                 fprintf(fp,"&\n");
 917             fprintf(fp,"%12.7f",time[i]);
 918             for(j=0; (j<nrms); j++)
 919                 fprintf(fp," %12.7f",rlsm[j][i]);
 920             fprintf(fp,"\n");
 921         }
 922         ffclose(fp);
 923     }
 924
 925     if (bAv) {
 926         sprintf(buf,"Average %s",whatxvgname[ewhat]);
 927         sprintf(buf2,"Average %s",whatxvglabel[ewhat]);
 928         fp = xvgropen(opt2fn("-a",NFILE,fnm), buf, "Residue", buf2,oenv);
 929         for(j=0; (j<nrms); j++)
 930             fprintf(fp,"%10d  %10g\n",j,rlstot/teller);
 931         ffclose(fp);
 932     }
 933
 934     if (bNorm) {
 935         fp = xvgropen("aver.xvg",gn_rms[0],"Residue",whatxvglabel[ewhat],oenv);
 936         for(j=0; (j<irms[0]); j++)
 937             fprintf(fp,"%10d  %10g\n",j,rlsnorm[j]/teller);
 938         ffclose(fp);
 939     }
 940     do_view(oenv,opt2fn_null("-a",NFILE,fnm),"-graphtype bar");
 941     do_view(oenv,opt2fn("-o",NFILE,fnm),NULL);
 942     do_view(oenv,opt2fn_null("-mir",NFILE,fnm),NULL);
 943     do_view(oenv,opt2fn_null("-m",NFILE,fnm),NULL);
 944     do_view(oenv,opt2fn_null("-bm",NFILE,fnm),NULL);
 945     do_view(oenv,opt2fn_null("-dist",NFILE,fnm),NULL);
 946
 947     thanx(stderr);
 948
 949     return 0;
 950 }