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     t_trxstatus *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     gmx_rmpbc_t  gpbc=NULL;
 249
 250     t_filenm fnm[] =
 251         {
 252             { efTPS, NULL, NULL, ffREAD },
 253             { efTRX, "-f", NULL, ffREAD },
 254             { efTRX, "-f2", NULL, ffOPTRD },
 255             { efNDX, NULL, NULL, ffOPTRD },
 256             { efXVG, NULL, "rmsd", ffWRITE },
 257             { efXVG, "-mir", "rmsdmir", ffOPTWR },
 258             { efXVG, "-a", "avgrp", ffOPTWR },
 259             { efXVG, "-dist", "rmsd-dist", ffOPTWR },
 260             { efXPM, "-m", "rmsd", ffOPTWR },
 261             { efDAT, "-bin", "rmsd", ffOPTWR },
 262             { efXPM, "-bm", "bond", ffOPTWR } };
 263 #define NFILE asize(fnm)
 264
 265     CopyRight(stderr, argv[0]);
 266     parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW
 267         | PCA_BE_NICE, NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL,
 268                       &oenv);
 269     /* parse enumerated options: */
 270     ewhat = nenum(what);
 271     if (ewhat == ewRho || ewhat == ewRhoSc)
 272         please_cite(stdout, "Maiorov95");
 273     efit = nenum(fit);
 274     bFit = efit == efFit;
 275     bReset = efit == efReset;
 276     if (bFit)
 277     {
 278         bReset = TRUE; /* for fit, reset *must* be set */
 279     }
 280     else
 281     {
 282         bFitAll = FALSE;
 283     }
 284
 285     /* mark active cmdline options */
 286     bMirror = opt2bSet("-mir", NFILE, fnm); /* calc RMSD vs mirror of ref. */
 287     bFile2 = opt2bSet("-f2", NFILE, fnm);
 288     bMat = opt2bSet("-m", NFILE, fnm);
 289     bBond = opt2bSet("-bm", NFILE, fnm);
 290     bDelta = (delta_maxy > 0); /* calculate rmsd vs delta t matrix from *
 291      *  your RMSD matrix (hidden option       */
 292     bNorm = opt2bSet("-a", NFILE, fnm);
 293     bFreq2 = opt2parg_bSet("-skip2", asize(pa), pa);
 294     if (freq <= 0)
 295     {
 296         fprintf(stderr, "The number of frames to skip is <= 0. "
 297             "Writing out all frames.\n\n");
 298         freq = 1;
 299     }
 300     if (!bFreq2)
 301     {
 302         freq2 = freq;
 303     }
 304     else if (bFile2 && freq2 <= 0)
 305     {
 306         fprintf(stderr,
 307                 "The number of frames to skip in second trajectory is <= 0.\n"
 308                     "  Writing out all frames.\n\n");
 309         freq2 = 1;
 310     }
 311
 312     bPrev = (prev > 0);
 313     if (bPrev)
 314     {
 315         prev = abs(prev);
 316         if (freq != 1)
 317             fprintf(stderr, "WARNING: option -skip also applies to -prev\n");
 318     }
 319
 320     if (bFile2 && !bMat && !bBond)
 321     {
 322         fprintf(
 323                 stderr,
 324                 "WARNING: second trajectory (-f2) useless when not calculating matrix (-m/-bm),\n"
 325                     "         will not read from %s\n", opt2fn("-f2", NFILE,
 326                                                                fnm));
 327         bFile2 = FALSE;
 328     }
 329
 330     if (bDelta)
 331     {
 332         bMat = TRUE;
 333         if (bFile2)
 334         {
 335             fprintf(stderr,
 336                     "WARNING: second trajectory (-f2) useless when making delta matrix,\n"
 337                         "         will not read from %s\n", opt2fn("-f2",
 338                                                                    NFILE, fnm));
 339             bFile2 = FALSE;
 340         }
 341     }
 342
 343     bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), buf, &top, &ePBC, &xp,
 344                          NULL, box, TRUE);
 345     snew(w_rls,top.atoms.nr);
 346     snew(w_rms,top.atoms.nr);
 347
 348     if (!bTop && bBond)
 349     {
 350         fprintf(stderr,
 351                 "WARNING: Need a run input file for bond angle matrix,\n"
 352                     "         will not calculate bond angle matrix.\n");
 353         bBond = FALSE;
 354     }
 355
 356     if (bReset)
 357     {
 358         fprintf(stderr, "Select group for %s fit\n", bFit ? "least squares"
 359             : "translational");
 360         get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE, fnm), 1, &ifit,
 361                   &ind_fit, &gn_fit);
 362     }
 363     else
 364         ifit = 0;
 365
 366     if (bReset)
 367     {
 368         if (bFit && ifit < 3)
 369             gmx_fatal(FARGS,"Need >= 3 points to fit!\n" );
 370
 371         bMass = FALSE;
 372         for(i=0; i<ifit; i++)
 373         {
 374             if (bMassWeighted)
 375                 w_rls[ind_fit[i]] = top.atoms.atom[ind_fit[i]].m;
 376             else
 377                 w_rls[ind_fit[i]] = 1;
 378             bMass = bMass || (top.atoms.atom[ind_fit[i]].m != 0);
 379         }
 380         if (!bMass)
 381         {
 382             fprintf(stderr,"All masses in the fit group are 0, using masses of 1\n");
 383             for(i=0; i<ifit; i++)
 384             {
 385                 w_rls[ind_fit[i]] = 1;
 386             }
 387         }
 388     }
 389
 390     if (bMat || bBond)
 391         nrms=1;
 392
 393     snew(gn_rms,nrms);
 394     snew(ind_rms,nrms);
 395     snew(irms,nrms);
 396
 397     fprintf(stderr,"Select group%s for %s calculation\n",
 398             (nrms>1) ? "s" : "",whatname[ewhat]);
 399     get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
 400               nrms,irms,ind_rms,gn_rms);
 401
 402     if (bNorm)
 403     {
 404         snew(rlsnorm,irms[0]);
 405     }
 406     snew(rls,nrms);
 407     for(j=0; j<nrms; j++)
 408         snew(rls[j],maxframe);
 409     if (bMirror)
 410     {
 411         snew(rlsm,nrms);
 412         for(j=0; j<nrms; j++)
 413             snew(rlsm[j],maxframe);
 414     }
 415     snew(time,maxframe);
 416     for(j=0; j<nrms; j++)
 417     {
 418         bMass = FALSE;
 419         for(i=0; i<irms[j]; i++)
 420         {
 421             if (bMassWeighted)
 422                 w_rms[ind_rms[j][i]] = top.atoms.atom[ind_rms[j][i]].m;
 423             else
 424                 w_rms[ind_rms[j][i]] = 1.0;
 425             bMass = bMass || (top.atoms.atom[ind_rms[j][i]].m != 0);
 426         }
 427         if (!bMass) {
 428             fprintf(stderr,"All masses in group %d are 0, using masses of 1\n",j);
 429             for(i=0; i<irms[j]; i++)
 430                 w_rms[ind_rms[j][i]] = 1;
 431         }
 432     }
 433     /* Prepare reference frame */
 434     if (bPBC) {
 435       gpbc = gmx_rmpbc_init(&top.idef,ePBC,top.atoms.nr,box);
 436       gmx_rmpbc(gpbc,top.atoms.nr,box,xp);
 437     }
 438     if (bReset)
 439         reset_x(ifit,ind_fit,top.atoms.nr,NULL,xp,w_rls);
 440     if (bMirror) {
 441         /* generate reference structure mirror image: */
 442         snew(xm, top.atoms.nr);
 443         for(i=0; i<top.atoms.nr; i++) {
 444             copy_rvec(xp[i],xm[i]);
 445             xm[i][XX] = -xm[i][XX];
 446         }
 447     }
 448     if (ewhat==ewRhoSc)
 449         norm_princ(&top.atoms, ifit, ind_fit, top.atoms.nr, xp);
 450
 451     /* read first frame */
 452     natoms=read_first_x(oenv,&status,opt2fn("-f",NFILE,fnm),&t,&x,box);
 453     if (natoms != top.atoms.nr)
 454         fprintf(stderr,
 455                 "\nWARNING: topology has %d atoms, whereas trajectory has %d\n",
 456                 top.atoms.nr,natoms);
 457     if (bMat || bBond || bPrev) {
 458         snew(mat_x,NFRAME);
 459
 460         if (bPrev)
 461             /* With -prev we use all atoms for simplicity */
 462             n_ind_m = natoms;
 463         else {
 464             /* Check which atoms we need (fit/rms) */
 465             snew(bInMat,natoms);
 466             for(i=0; i<ifit; i++)
 467                 bInMat[ind_fit[i]] = TRUE;
 468             n_ind_m=ifit;
 469             for(i=0; i<irms[0]; i++)
 470                 if (!bInMat[ind_rms[0][i]]) {
 471                     bInMat[ind_rms[0][i]] = TRUE;
 472                     n_ind_m++;
 473                 }
 474         }
 475         /* Make an index of needed atoms */
 476         snew(ind_m,n_ind_m);
 477         snew(rev_ind_m,natoms);
 478         j = 0;
 479         for(i=0; i<natoms; i++)
 480             if (bPrev || bInMat[i]) {
 481                 ind_m[j] = i;
 482                 rev_ind_m[i] = j;
 483                 j++;
 484             }
 485         snew(w_rls_m,n_ind_m);
 486         snew(ind_rms_m,irms[0]);
 487         snew(w_rms_m,n_ind_m);
 488         for(i=0; i<ifit; i++)
 489             w_rls_m[rev_ind_m[ind_fit[i]]] = w_rls[ind_fit[i]];
 490         for(i=0; i<irms[0]; i++) {
 491             ind_rms_m[i] = rev_ind_m[ind_rms[0][i]];
 492             w_rms_m[ind_rms_m[i]] = w_rms[ind_rms[0][i]];
 493         }
 494         sfree(bInMat);
 495     }
 496
 497     if (bBond) {
 498         ncons = 0;
 499         for(k=0; k<F_NRE; k++)
 500             if (IS_CHEMBOND(k)) {
 501                 iatom  = top.idef.il[k].iatoms;
 502                 ncons += top.idef.il[k].nr/3;
 503             }
 504         fprintf(stderr,"Found %d bonds in topology\n",ncons);
 505         snew(ind_bond1,ncons);
 506         snew(ind_bond2,ncons);
 507         ibond=0;
 508         for(k=0; k<F_NRE; k++)
 509             if (IS_CHEMBOND(k)) {
 510                 iatom = top.idef.il[k].iatoms;
 511                 ncons = top.idef.il[k].nr/3;
 512                 for (i=0; i<ncons; i++) {
 513                     bA1=FALSE;
 514                     bA2=FALSE;
 515                     for (j=0; j<irms[0]; j++) {
 516                         if (iatom[3*i+1] == ind_rms[0][j]) bA1=TRUE;
 517                         if (iatom[3*i+2] == ind_rms[0][j]) bA2=TRUE;
 518                     }
 519                     if (bA1 && bA2) {
 520                         ind_bond1[ibond] = rev_ind_m[iatom[3*i+1]];
 521                         ind_bond2[ibond] = rev_ind_m[iatom[3*i+2]];
 522                         ibond++;
 523                     }
 524                 }
 525             }
 526         fprintf(stderr,"Using %d bonds for bond angle matrix\n",ibond);
 527         if (ibond==0)
 528             gmx_fatal(FARGS,"0 bonds found");
 529     }
 530
 531     /* start looping over frames: */
 532     tel_mat = 0;
 533     teller = 0;
 534     do {
 535         if (bPBC)
 536           gmx_rmpbc(gpbc,natoms,box,x);
 537
 538         if (bReset)
 539             reset_x(ifit,ind_fit,natoms,NULL,x,w_rls);
 540         if (ewhat==ewRhoSc)
 541             norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 542
 543         if (bFit)
 544             /*do the least squares fit to original structure*/
 545             do_fit(natoms,w_rls,xp,x);
 546
 547         if (teller % freq == 0) {
 548             /* keep frame for matrix calculation */
 549             if (bMat || bBond || bPrev) {
 550                 if (tel_mat >= NFRAME)
 551                     srenew(mat_x,tel_mat+1);
 552                 snew(mat_x[tel_mat],n_ind_m);
 553                 for (i=0;i<n_ind_m;i++)
 554                     copy_rvec(x[ind_m[i]],mat_x[tel_mat][i]);
 555             }
 556             tel_mat++;
 557         }
 558
 559         /*calculate energy of root_least_squares*/
 560         if (bPrev) {
 561             j=tel_mat-prev-1;
 562             if (j<0)
 563                 j=0;
 564             for (i=0;i<n_ind_m;i++)
 565                 copy_rvec(mat_x[j][i],xp[ind_m[i]]);
 566             if (bReset)
 567                 reset_x(ifit,ind_fit,natoms,NULL,xp,w_rls);
 568             if (bFit)
 569                 do_fit(natoms,w_rls,x,xp);
 570         }
 571         for(j=0; (j<nrms); j++)
 572             rls[j][teller] =
 573                 calc_similar_ind(ewhat!=ewRMSD,irms[j],ind_rms[j],w_rms,x,xp);
 574         if (bNorm) {
 575             for(j=0; (j<irms[0]); j++)
 576                 rlsnorm[j] +=
 577                     calc_similar_ind(ewhat!=ewRMSD,1,&(ind_rms[0][j]),w_rms,x,xp);
 578         }
 579
 580         if (bMirror) {
 581             if (bFit)
 582                 /*do the least squares fit to mirror of original structure*/
 583                 do_fit(natoms,w_rls,xm,x);
 584
 585             for(j=0; j<nrms; j++)
 586                 rlsm[j][teller] =
 587                     calc_similar_ind(ewhat!=ewRMSD,irms[j],ind_rms[j],w_rms,x,xm);
 588         }
 589         time[teller]=output_env_conv_time(oenv,t);
 590
 591         teller++;
 592         if (teller >= maxframe) {
 593             maxframe +=NFRAME;
 594             srenew(time,maxframe);
 595             for(j=0; (j<nrms); j++)
 596                 srenew(rls[j],maxframe);
 597             if (bMirror)
 598                 for(j=0; (j<nrms); j++)
 599                     srenew(rlsm[j],maxframe);
 600         }
 601     } while (read_next_x(oenv,status,&t,natoms,x,box));
 602     close_trj(status);
 603
 604     if (bFile2) {
 605         snew(time2,maxframe2);
 606
 607         fprintf(stderr,"\nWill read second trajectory file\n");
 608         snew(mat_x2,NFRAME);
 609         natoms2=read_first_x(oenv,&status,opt2fn("-f2",NFILE,fnm),&t,&x,box);
 610         if ( natoms2 != natoms )
 611             gmx_fatal(FARGS,
 612                       "Second trajectory (%d atoms) does not match the first one"
 613                       " (%d atoms)", natoms2, natoms);
 614         tel_mat2 = 0;
 615         teller2 = 0;
 616         do {
 617             if (bPBC)
 618               gmx_rmpbc(gpbc,natoms,box,x);
 619
 620             if (bReset)
 621                 reset_x(ifit,ind_fit,natoms,NULL,x,w_rls);
 622             if (ewhat==ewRhoSc)
 623                 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 624
 625             if (bFit)
 626                 /*do the least squares fit to original structure*/
 627                 do_fit(natoms,w_rls,xp,x);
 628
 629             if (teller2 % freq2 == 0) {
 630                 /* keep frame for matrix calculation */
 631                 if (bMat) {
 632                     if (tel_mat2 >= NFRAME)
 633                         srenew(mat_x2,tel_mat2+1);
 634                     snew(mat_x2[tel_mat2],n_ind_m);
 635                     for (i=0;i<n_ind_m;i++)
 636                         copy_rvec(x[ind_m[i]],mat_x2[tel_mat2][i]);
 637                 }
 638                 tel_mat2++;
 639             }
 640
 641             time2[teller2]=output_env_conv_time(oenv,t);
 642
 643             teller2++;
 644             if (teller2 >= maxframe2) {
 645                 maxframe2 +=NFRAME;
 646                 srenew(time2,maxframe2);
 647             }
 648         } while (read_next_x(oenv,status,&t,natoms,x,box));
 649         close_trj(status);
 650     } else {
 651         mat_x2=mat_x;
 652         time2=time;
 653         tel_mat2=tel_mat;
 654         freq2=freq;
 655     }
 656     gmx_rmpbc_done(gpbc);
 657
 658     if (bMat || bBond) {
 659         /* calculate RMS matrix */
 660         fprintf(stderr,"\n");
 661         if (bMat) {
 662             fprintf(stderr,"Building %s matrix, %dx%d elements\n",
 663                     whatname[ewhat],tel_mat,tel_mat2);
 664             snew(rmsd_mat,tel_mat);
 665         }
 666         if (bBond) {
 667             fprintf(stderr,"Building bond angle matrix, %dx%d elements\n",
 668                     tel_mat,tel_mat2);
 669             snew(bond_mat,tel_mat);
 670         }
 671         snew(axis,tel_mat);
 672         snew(axis2,tel_mat2);
 673         rmsd_max=0;
 674         if (bFile2)
 675             rmsd_min=1e10;
 676         else
 677             rmsd_min=0;
 678         rmsd_avg=0;
 679         bond_max=0;
 680         bond_min=1e10;
 681         for(j=0; j<tel_mat2; j++)
 682             axis2[j]=time2[freq2*j];
 683         if (bDelta) {
 684             if (bDeltaLog) {
 685                 delta_scalex=8.0/log(2.0);
 686                 delta_xsize=(int)(log(tel_mat/2)*delta_scalex+0.5)+1;
 687             }
 688             else {
 689                 delta_xsize=tel_mat/2;
 690             }
 691             delta_scaley=1.0/delta_maxy;
 692             snew(delta,delta_xsize);
 693             for(j=0; j<delta_xsize; j++)
 694                 snew(delta[j],del_lev+1);
 695             if (avl > 0) {
 696                 snew(rmsdav_mat,tel_mat);
 697                 for(j=0; j<tel_mat; j++)
 698                     snew(rmsdav_mat[j],tel_mat);
 699             }
 700         }
 701
 702         if (bFitAll)
 703             snew(mat_x2_j,natoms);
 704         for(i=0; i<tel_mat; i++) {
 705             axis[i]=time[freq*i];
 706             fprintf(stderr,"\r element %5d; time %5.2f  ",i,axis[i]);
 707             if (bMat) snew(rmsd_mat[i],tel_mat2);
 708             if (bBond) snew(bond_mat[i],tel_mat2);
 709             for(j=0; j<tel_mat2; j++) {
 710                 if (bFitAll) {
 711                     for (k=0;k<n_ind_m;k++)
 712                         copy_rvec(mat_x2[j][k],mat_x2_j[k]);
 713                     do_fit(n_ind_m,w_rls_m,mat_x[i],mat_x2_j);
 714                 } else
 715                     mat_x2_j=mat_x2[j];
 716                 if (bMat) {
 717                     if (bFile2 || (i<j)) {
 718                         rmsd_mat[i][j] =
 719                             calc_similar_ind(ewhat!=ewRMSD,irms[0],ind_rms_m,
 720                                              w_rms_m,mat_x[i],mat_x2_j);
 721                         if (rmsd_mat[i][j] > rmsd_max) rmsd_max=rmsd_mat[i][j];
 722                         if (rmsd_mat[i][j] < rmsd_min) rmsd_min=rmsd_mat[i][j];
 723                         rmsd_avg += rmsd_mat[i][j];
 724                     } else
 725                         rmsd_mat[i][j]=rmsd_mat[j][i];
 726                 }
 727                 if (bBond) {
 728                     if (bFile2 || (i<=j)) {
 729                         ang=0.0;
 730                         for(m=0;m<ibond;m++) {
 731                             rvec_sub(mat_x[i][ind_bond1[m]],mat_x[i][ind_bond2[m]],vec1);
 732                             rvec_sub(mat_x2_j[ind_bond1[m]],mat_x2_j[ind_bond2[m]],vec2);
 733                             ang += acos(cos_angle(vec1,vec2));
 734                         }
 735                         bond_mat[i][j]=ang*180.0/(M_PI*ibond);
 736                         if (bond_mat[i][j] > bond_max) bond_max=bond_mat[i][j];
 737                         if (bond_mat[i][j] < bond_min) bond_min=bond_mat[i][j];
 738                     }
 739                     else
 740                         bond_mat[i][j]=bond_mat[j][i];
 741                 }
 742             }
 743         }
 744         if (bFile2)
 745             rmsd_avg /= tel_mat*tel_mat2;
 746         else
 747             rmsd_avg /= tel_mat*(tel_mat - 1)/2;
 748         if (bMat && (avl > 0)) {
 749             rmsd_max=0.0;
 750             rmsd_min=0.0;
 751             rmsd_avg=0.0;
 752             for(j=0; j<tel_mat-1; j++) {
 753                 for(i=j+1; i<tel_mat; i++) {
 754                     av_tot=0;
 755                     weight_tot=0;
 756                     for (my=-avl; my<=avl; my++) {
 757                         if ((j+my>=0) && (j+my<tel_mat)) {
 758                             abs_my = abs(my);
 759                             for (mx=-avl; mx<=avl; mx++) {
 760                                 if ((i+mx>=0) && (i+mx<tel_mat)) {
 761                                     weight = (real)(avl+1-max(abs(mx),abs_my));
 762                                     av_tot += weight*rmsd_mat[i+mx][j+my];
 763                                     weight_tot+=weight;
 764                                 }
 765                             }
 766                         }
 767                     }
 768                     rmsdav_mat[i][j] = av_tot/weight_tot;
 769                     rmsdav_mat[j][i] = rmsdav_mat[i][j];
 770                     if (rmsdav_mat[i][j] > rmsd_max) rmsd_max=rmsdav_mat[i][j];
 771                 }
 772             }
 773             rmsd_mat=rmsdav_mat;
 774         }
 775
 776         if (bMat) {
 777             fprintf(stderr,"\n%s: Min %f, Max %f, Avg %f\n",
 778                     whatname[ewhat],rmsd_min,rmsd_max,rmsd_avg);
 779             rlo.r = 1; rlo.g = 1; rlo.b = 1;
 780             rhi.r = 0; rhi.g = 0; rhi.b = 0;
 781             if (rmsd_user_max != -1) rmsd_max=rmsd_user_max;
 782             if (rmsd_user_min != -1) rmsd_min=rmsd_user_min;
 783             if ((rmsd_user_max !=  -1) || (rmsd_user_min != -1))
 784                 fprintf(stderr,"Min and Max value set to resp. %f and %f\n",
 785                         rmsd_min,rmsd_max);
 786             sprintf(buf,"%s %s matrix",gn_rms[0],whatname[ewhat]);
 787             write_xpm(opt2FILE("-m",NFILE,fnm,"w"),0,buf,whatlabel[ewhat],
 788                       output_env_get_time_label(oenv),output_env_get_time_label(oenv),tel_mat,tel_mat2,
 789                       axis,axis2, rmsd_mat,rmsd_min,rmsd_max,rlo,rhi,&nlevels);
 790             /* Print the distribution of RMSD values */
 791             if (opt2bSet("-dist",NFILE,fnm))
 792                 low_rmsd_dist(opt2fn("-dist",NFILE,fnm),rmsd_max,tel_mat,rmsd_mat,oenv);
 793
 794             if (bDelta) {
 795                 snew(delta_tot,delta_xsize);
 796                 for(j=0; j<tel_mat-1; j++) {
 797                     for(i=j+1; i<tel_mat; i++) {
 798                         mx=i-j ;
 799                         if (mx < tel_mat/2) {
 800                             if (bDeltaLog)
 801                                 mx=(int)(log(mx)*delta_scalex+0.5);
 802                             my=(int)(rmsd_mat[i][j]*delta_scaley*del_lev+0.5);
 803                             delta_tot[mx] += 1.0;
 804                             if ((rmsd_mat[i][j]>=0) && (rmsd_mat[i][j]<=delta_maxy))
 805                                 delta[mx][my] += 1.0;
 806                         }
 807                     }
 808                 }
 809                 delta_max=0;
 810                 for(i=0; i<delta_xsize; i++) {
 811                     if (delta_tot[i] > 0.0) {
 812                         delta_tot[i]=1.0/delta_tot[i];
 813                         for(j=0; j<=del_lev; j++) {
 814                             delta[i][j] *= delta_tot[i];
 815                             if (delta[i][j] > delta_max)
 816                                 delta_max=delta[i][j];
 817                         }
 818                     }
 819                 }
 820                 fprintf(stderr,"Maximum in delta matrix: %f\n",delta_max);
 821                 snew(del_xaxis,delta_xsize);
 822                 snew(del_yaxis,del_lev+1);
 823                 for (i=0; i<delta_xsize; i++)
 824                     del_xaxis[i]=axis[i]-axis[0];
 825                 for (i=0; i<del_lev+1; i++)
 826                     del_yaxis[i]=delta_maxy*i/del_lev;
 827                 sprintf(buf,"%s %s vs. delta t",gn_rms[0],whatname[ewhat]);
 828                 fp = ffopen("delta.xpm","w");
 829                 write_xpm(fp,0,buf,"density",output_env_get_time_label(oenv),whatlabel[ewhat],
 830                           delta_xsize,del_lev+1,del_xaxis,del_yaxis,
 831                           delta,0.0,delta_max,rlo,rhi,&nlevels);
 832                 ffclose(fp);
 833             }
 834             if (opt2bSet("-bin",NFILE,fnm)) {
 835                 /* NB: File must be binary if we use fwrite */
 836                 fp=ftp2FILE(efDAT,NFILE,fnm,"wb");
 837                 for(i=0;i<tel_mat;i++)
 838                     if(fwrite(rmsd_mat[i],sizeof(**rmsd_mat),tel_mat2,fp) != tel_mat2)
 839                     {
 840                         gmx_fatal(FARGS,"Error writing to output file");
 841                     }
 842                 ffclose(fp);
 843             }
 844         }
 845         if (bBond) {
 846             fprintf(stderr,"\nMin. angle: %f, Max. angle: %f\n",bond_min,bond_max);
 847             if (bond_user_max != -1) bond_max=bond_user_max;
 848             if (bond_user_min != -1) bond_min=bond_user_min;
 849             if ((bond_user_max !=  -1) || (bond_user_min != -1))
 850                 fprintf(stderr,"Bond angle Min and Max set to:\n"
 851                         "Min. angle: %f, Max. angle: %f\n",bond_min,bond_max);
 852             rlo.r = 1; rlo.g = 1; rlo.b = 1;
 853             rhi.r = 0; rhi.g = 0; rhi.b = 0;
 854             sprintf(buf,"%s av. bond angle deviation",gn_rms[0]);
 855             write_xpm(opt2FILE("-bm",NFILE,fnm,"w"),0,buf,"degrees",
 856                       output_env_get_time_label(oenv),output_env_get_time_label(oenv),tel_mat,tel_mat2,
 857                       axis,axis2, bond_mat,bond_min,bond_max,rlo,rhi,&nlevels);
 858         }
 859     }
 860
 861     bAv=opt2bSet("-a",NFILE,fnm);
 862
 863     /* Write the RMSD's to file */
 864     if (!bPrev)
 865         sprintf(buf,"%s",whatxvgname[ewhat]);
 866     else
 867         sprintf(buf,"%s with frame %g %s ago",whatxvgname[ewhat],
 868                 time[prev*freq]-time[0], output_env_get_time_label(oenv));
 869     fp=xvgropen(opt2fn("-o",NFILE,fnm),buf,output_env_get_xvgr_tlabel(oenv),
 870                 whatxvglabel[ewhat], oenv);
 871     if (output_env_get_print_xvgr_codes(oenv))
 872         fprintf(fp,"@ subtitle \"%s%s after %s%s%s\"\n",
 873                 (nrms==1)?"":"of "    , gn_rms[0], fitgraphlabel[efit],
 874                     bFit     ?" to ":""   , bFit?gn_fit:"");
 875     if (nrms != 1)
 876         xvgr_legend(fp,nrms,(const char**)gn_rms,oenv);
 877     for(i=0; (i<teller); i++) {
 878         if ( bSplit && i>0 &&
 879             abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv))<1e-5 )
 880         {
 881             fprintf(fp,"&\n");
 882         }
 883         fprintf(fp,"%12.7f",time[bPrev ? freq*i : i]);
 884         for(j=0; (j<nrms); j++) {
 885             fprintf(fp," %12.7f",rls[j][i]);
 886             if (bAv)
 887                 rlstot+=rls[j][i];
 888         }
 889         fprintf(fp,"\n");
 890     }
 891     ffclose(fp);
 892
 893     if (bMirror) {
 894         /* Write the mirror RMSD's to file */
 895         sprintf(buf,"%s with Mirror",whatxvgname[ewhat]);
 896         sprintf(buf2,"Mirror %s",whatxvglabel[ewhat]);
 897         fp=xvgropen(opt2fn("-mir",NFILE,fnm), buf, output_env_get_xvgr_tlabel(oenv),
 898                     buf2,oenv);
 899         if (nrms == 1) {
 900             if (output_env_get_print_xvgr_codes(oenv))
 901                 fprintf(fp,"@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
 902                         gn_rms[0],gn_fit);
 903         }
 904         else {
 905             if (output_env_get_print_xvgr_codes(oenv))
 906                 fprintf(fp,"@ subtitle \"after lsq fit to mirror %s\"\n",gn_fit);
 907             xvgr_legend(fp,nrms,(const char**)gn_rms,oenv);
 908         }
 909         for(i=0; (i<teller); i++) {
 910             if ( bSplit && i>0 && abs(time[i])<1e-5 )
 911                 fprintf(fp,"&\n");
 912             fprintf(fp,"%12.7f",time[i]);
 913             for(j=0; (j<nrms); j++)
 914                 fprintf(fp," %12.7f",rlsm[j][i]);
 915             fprintf(fp,"\n");
 916         }
 917         ffclose(fp);
 918     }
 919
 920     if (bAv) {
 921         sprintf(buf,"Average %s",whatxvgname[ewhat]);
 922         sprintf(buf2,"Average %s",whatxvglabel[ewhat]);
 923         fp = xvgropen(opt2fn("-a",NFILE,fnm), buf, "Residue", buf2,oenv);
 924         for(j=0; (j<nrms); j++)
 925             fprintf(fp,"%10d  %10g\n",j,rlstot/teller);
 926         ffclose(fp);
 927     }
 928
 929     if (bNorm) {
 930         fp = xvgropen("aver.xvg",gn_rms[0],"Residue",whatxvglabel[ewhat],oenv);
 931         for(j=0; (j<irms[0]); j++)
 932             fprintf(fp,"%10d  %10g\n",j,rlsnorm[j]/teller);
 933         ffclose(fp);
 934     }
 935     do_view(oenv,opt2fn_null("-a",NFILE,fnm),"-graphtype bar");
 936     do_view(oenv,opt2fn("-o",NFILE,fnm),NULL);
 937     do_view(oenv,opt2fn_null("-mir",NFILE,fnm),NULL);
 938     do_view(oenv,opt2fn_null("-m",NFILE,fnm),NULL);
 939     do_view(oenv,opt2fn_null("-bm",NFILE,fnm),NULL);
 940     do_view(oenv,opt2fn_null("-dist",NFILE,fnm),NULL);
 941
 942     thanx(stderr);
 943
 944     return 0;
 945 }