src/tools/gmx_rms.c

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