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                 "[TT]g_rms[tt] 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 [TT].tpr[tt] file masses will be taken from there, ",
 127                 "otherwise the masses will be deduced from the [TT]atommass.dat[tt] file in",
 128                 "[TT]GMXLIB[tt]. 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 gmx_bool bPBC = TRUE, bFitAll = TRUE, bSplit = FALSE;
 143     static gmx_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 gmx_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_trx, natoms_trx2, natoms;
 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     gmx_bool bNorm, bAv, bFreq2, bFile2, bMat, bBond, bDelta, bMirror, bMass;
 221     gmx_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     gmx_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_trx=read_first_x(oenv,&status,opt2fn("-f",NFILE,fnm),&t,&x,box);
 453     if (natoms_trx != top.atoms.nr)
 454         fprintf(stderr,
 455                 "\nWARNING: topology has %d atoms, whereas trajectory has %d\n",
 456                 top.atoms.nr,natoms_trx);
 457     natoms = min(top.atoms.nr,natoms_trx);
 458     if (bMat || bBond || bPrev) {
 459         snew(mat_x,NFRAME);
 460
 461         if (bPrev)
 462             /* With -prev we use all atoms for simplicity */
 463             n_ind_m = natoms;
 464         else {
 465             /* Check which atoms we need (fit/rms) */
 466             snew(bInMat,natoms);
 467             for(i=0; i<ifit; i++)
 468                 bInMat[ind_fit[i]] = TRUE;
 469             n_ind_m=ifit;
 470             for(i=0; i<irms[0]; i++)
 471                 if (!bInMat[ind_rms[0][i]]) {
 472                     bInMat[ind_rms[0][i]] = TRUE;
 473                     n_ind_m++;
 474                 }
 475         }
 476         /* Make an index of needed atoms */
 477         snew(ind_m,n_ind_m);
 478         snew(rev_ind_m,natoms);
 479         j = 0;
 480         for(i=0; i<natoms; i++)
 481             if (bPrev || bInMat[i]) {
 482                 ind_m[j] = i;
 483                 rev_ind_m[i] = j;
 484                 j++;
 485             }
 486         snew(w_rls_m,n_ind_m);
 487         snew(ind_rms_m,irms[0]);
 488         snew(w_rms_m,n_ind_m);
 489         for(i=0; i<ifit; i++)
 490             w_rls_m[rev_ind_m[ind_fit[i]]] = w_rls[ind_fit[i]];
 491         for(i=0; i<irms[0]; i++) {
 492             ind_rms_m[i] = rev_ind_m[ind_rms[0][i]];
 493             w_rms_m[ind_rms_m[i]] = w_rms[ind_rms[0][i]];
 494         }
 495         sfree(bInMat);
 496     }
 497
 498     if (bBond) {
 499         ncons = 0;
 500         for(k=0; k<F_NRE; k++)
 501             if (IS_CHEMBOND(k)) {
 502                 iatom  = top.idef.il[k].iatoms;
 503                 ncons += top.idef.il[k].nr/3;
 504             }
 505         fprintf(stderr,"Found %d bonds in topology\n",ncons);
 506         snew(ind_bond1,ncons);
 507         snew(ind_bond2,ncons);
 508         ibond=0;
 509         for(k=0; k<F_NRE; k++)
 510             if (IS_CHEMBOND(k)) {
 511                 iatom = top.idef.il[k].iatoms;
 512                 ncons = top.idef.il[k].nr/3;
 513                 for (i=0; i<ncons; i++) {
 514                     bA1=FALSE;
 515                     bA2=FALSE;
 516                     for (j=0; j<irms[0]; j++) {
 517                         if (iatom[3*i+1] == ind_rms[0][j]) bA1=TRUE;
 518                         if (iatom[3*i+2] == ind_rms[0][j]) bA2=TRUE;
 519                     }
 520                     if (bA1 && bA2) {
 521                         ind_bond1[ibond] = rev_ind_m[iatom[3*i+1]];
 522                         ind_bond2[ibond] = rev_ind_m[iatom[3*i+2]];
 523                         ibond++;
 524                     }
 525                 }
 526             }
 527         fprintf(stderr,"Using %d bonds for bond angle matrix\n",ibond);
 528         if (ibond==0)
 529             gmx_fatal(FARGS,"0 bonds found");
 530     }
 531
 532     /* start looping over frames: */
 533     tel_mat = 0;
 534     teller = 0;
 535     do {
 536         if (bPBC)
 537           gmx_rmpbc(gpbc,natoms,box,x);
 538
 539         if (bReset)
 540             reset_x(ifit,ind_fit,natoms,NULL,x,w_rls);
 541         if (ewhat==ewRhoSc)
 542             norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 543
 544         if (bFit)
 545             /*do the least squares fit to original structure*/
 546             do_fit(natoms,w_rls,xp,x);
 547
 548         if (teller % freq == 0) {
 549             /* keep frame for matrix calculation */
 550             if (bMat || bBond || bPrev) {
 551                 if (tel_mat >= NFRAME)
 552                     srenew(mat_x,tel_mat+1);
 553                 snew(mat_x[tel_mat],n_ind_m);
 554                 for (i=0;i<n_ind_m;i++)
 555                     copy_rvec(x[ind_m[i]],mat_x[tel_mat][i]);
 556             }
 557             tel_mat++;
 558         }
 559
 560         /*calculate energy of root_least_squares*/
 561         if (bPrev) {
 562             j=tel_mat-prev-1;
 563             if (j<0)
 564                 j=0;
 565             for (i=0;i<n_ind_m;i++)
 566                 copy_rvec(mat_x[j][i],xp[ind_m[i]]);
 567             if (bReset)
 568                 reset_x(ifit,ind_fit,natoms,NULL,xp,w_rls);
 569             if (bFit)
 570                 do_fit(natoms,w_rls,x,xp);
 571         }
 572         for(j=0; (j<nrms); j++)
 573             rls[j][teller] =
 574                 calc_similar_ind(ewhat!=ewRMSD,irms[j],ind_rms[j],w_rms,x,xp);
 575         if (bNorm) {
 576             for(j=0; (j<irms[0]); j++)
 577                 rlsnorm[j] +=
 578                     calc_similar_ind(ewhat!=ewRMSD,1,&(ind_rms[0][j]),w_rms,x,xp);
 579         }
 580
 581         if (bMirror) {
 582             if (bFit)
 583                 /*do the least squares fit to mirror of original structure*/
 584                 do_fit(natoms,w_rls,xm,x);
 585
 586             for(j=0; j<nrms; j++)
 587                 rlsm[j][teller] =
 588                     calc_similar_ind(ewhat!=ewRMSD,irms[j],ind_rms[j],w_rms,x,xm);
 589         }
 590         time[teller]=output_env_conv_time(oenv,t);
 591
 592         teller++;
 593         if (teller >= maxframe) {
 594             maxframe +=NFRAME;
 595             srenew(time,maxframe);
 596             for(j=0; (j<nrms); j++)
 597                 srenew(rls[j],maxframe);
 598             if (bMirror)
 599                 for(j=0; (j<nrms); j++)
 600                     srenew(rlsm[j],maxframe);
 601         }
 602     } while (read_next_x(oenv,status,&t,natoms_trx,x,box));
 603     close_trj(status);
 604
 605     if (bFile2) {
 606         snew(time2,maxframe2);
 607
 608         fprintf(stderr,"\nWill read second trajectory file\n");
 609         snew(mat_x2,NFRAME);
 610         natoms_trx2 =
 611           read_first_x(oenv,&status,opt2fn("-f2",NFILE,fnm),&t,&x,box);
 612         if ( natoms_trx2 != natoms_trx )
 613             gmx_fatal(FARGS,
 614                       "Second trajectory (%d atoms) does not match the first one"
 615                       " (%d atoms)", natoms_trx2, natoms_trx);
 616         tel_mat2 = 0;
 617         teller2 = 0;
 618         do {
 619             if (bPBC)
 620               gmx_rmpbc(gpbc,natoms,box,x);
 621
 622             if (bReset)
 623                 reset_x(ifit,ind_fit,natoms,NULL,x,w_rls);
 624             if (ewhat==ewRhoSc)
 625                 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
 626
 627             if (bFit)
 628                 /*do the least squares fit to original structure*/
 629                 do_fit(natoms,w_rls,xp,x);
 630
 631             if (teller2 % freq2 == 0) {
 632                 /* keep frame for matrix calculation */
 633                 if (bMat) {
 634                     if (tel_mat2 >= NFRAME)
 635                         srenew(mat_x2,tel_mat2+1);
 636                     snew(mat_x2[tel_mat2],n_ind_m);
 637                     for (i=0;i<n_ind_m;i++)
 638                         copy_rvec(x[ind_m[i]],mat_x2[tel_mat2][i]);
 639                 }
 640                 tel_mat2++;
 641             }
 642
 643             time2[teller2]=output_env_conv_time(oenv,t);
 644
 645             teller2++;
 646             if (teller2 >= maxframe2) {
 647                 maxframe2 +=NFRAME;
 648                 srenew(time2,maxframe2);
 649             }
 650         } while (read_next_x(oenv,status,&t,natoms_trx2,x,box));
 651         close_trj(status);
 652     } else {
 653         mat_x2=mat_x;
 654         time2=time;
 655         tel_mat2=tel_mat;
 656         freq2=freq;
 657     }
 658     gmx_rmpbc_done(gpbc);
 659
 660     if (bMat || bBond) {
 661         /* calculate RMS matrix */
 662         fprintf(stderr,"\n");
 663         if (bMat) {
 664             fprintf(stderr,"Building %s matrix, %dx%d elements\n",
 665                     whatname[ewhat],tel_mat,tel_mat2);
 666             snew(rmsd_mat,tel_mat);
 667         }
 668         if (bBond) {
 669             fprintf(stderr,"Building bond angle matrix, %dx%d elements\n",
 670                     tel_mat,tel_mat2);
 671             snew(bond_mat,tel_mat);
 672         }
 673         snew(axis,tel_mat);
 674         snew(axis2,tel_mat2);
 675         rmsd_max=0;
 676         if (bFile2)
 677             rmsd_min=1e10;
 678         else
 679             rmsd_min=0;
 680         rmsd_avg=0;
 681         bond_max=0;
 682         bond_min=1e10;
 683         for(j=0; j<tel_mat2; j++)
 684             axis2[j]=time2[freq2*j];
 685         if (bDelta) {
 686             if (bDeltaLog) {
 687                 delta_scalex=8.0/log(2.0);
 688                 delta_xsize=(int)(log(tel_mat/2)*delta_scalex+0.5)+1;
 689             }
 690             else {
 691                 delta_xsize=tel_mat/2;
 692             }
 693             delta_scaley=1.0/delta_maxy;
 694             snew(delta,delta_xsize);
 695             for(j=0; j<delta_xsize; j++)
 696                 snew(delta[j],del_lev+1);
 697             if (avl > 0) {
 698                 snew(rmsdav_mat,tel_mat);
 699                 for(j=0; j<tel_mat; j++)
 700                     snew(rmsdav_mat[j],tel_mat);
 701             }
 702         }
 703
 704         if (bFitAll)
 705             snew(mat_x2_j,natoms);
 706         for(i=0; i<tel_mat; i++) {
 707             axis[i]=time[freq*i];
 708             fprintf(stderr,"\r element %5d; time %5.2f  ",i,axis[i]);
 709             if (bMat) snew(rmsd_mat[i],tel_mat2);
 710             if (bBond) snew(bond_mat[i],tel_mat2);
 711             for(j=0; j<tel_mat2; j++) {
 712                 if (bFitAll) {
 713                     for (k=0;k<n_ind_m;k++)
 714                         copy_rvec(mat_x2[j][k],mat_x2_j[k]);
 715                     do_fit(n_ind_m,w_rls_m,mat_x[i],mat_x2_j);
 716                 } else
 717                     mat_x2_j=mat_x2[j];
 718                 if (bMat) {
 719                     if (bFile2 || (i<j)) {
 720                         rmsd_mat[i][j] =
 721                             calc_similar_ind(ewhat!=ewRMSD,irms[0],ind_rms_m,
 722                                              w_rms_m,mat_x[i],mat_x2_j);
 723                         if (rmsd_mat[i][j] > rmsd_max) rmsd_max=rmsd_mat[i][j];
 724                         if (rmsd_mat[i][j] < rmsd_min) rmsd_min=rmsd_mat[i][j];
 725                         rmsd_avg += rmsd_mat[i][j];
 726                     } else
 727                         rmsd_mat[i][j]=rmsd_mat[j][i];
 728                 }
 729                 if (bBond) {
 730                     if (bFile2 || (i<=j)) {
 731                         ang=0.0;
 732                         for(m=0;m<ibond;m++) {
 733                             rvec_sub(mat_x[i][ind_bond1[m]],mat_x[i][ind_bond2[m]],vec1);
 734                             rvec_sub(mat_x2_j[ind_bond1[m]],mat_x2_j[ind_bond2[m]],vec2);
 735                             ang += acos(cos_angle(vec1,vec2));
 736                         }
 737                         bond_mat[i][j]=ang*180.0/(M_PI*ibond);
 738                         if (bond_mat[i][j] > bond_max) bond_max=bond_mat[i][j];
 739                         if (bond_mat[i][j] < bond_min) bond_min=bond_mat[i][j];
 740                     }
 741                     else
 742                         bond_mat[i][j]=bond_mat[j][i];
 743                 }
 744             }
 745         }
 746         if (bFile2)
 747             rmsd_avg /= tel_mat*tel_mat2;
 748         else
 749             rmsd_avg /= tel_mat*(tel_mat - 1)/2;
 750         if (bMat && (avl > 0)) {
 751             rmsd_max=0.0;
 752             rmsd_min=0.0;
 753             rmsd_avg=0.0;
 754             for(j=0; j<tel_mat-1; j++) {
 755                 for(i=j+1; i<tel_mat; i++) {
 756                     av_tot=0;
 757                     weight_tot=0;
 758                     for (my=-avl; my<=avl; my++) {
 759                         if ((j+my>=0) && (j+my<tel_mat)) {
 760                             abs_my = abs(my);
 761                             for (mx=-avl; mx<=avl; mx++) {
 762                                 if ((i+mx>=0) && (i+mx<tel_mat)) {
 763                                     weight = (real)(avl+1-max(abs(mx),abs_my));
 764                                     av_tot += weight*rmsd_mat[i+mx][j+my];
 765                                     weight_tot+=weight;
 766                                 }
 767                             }
 768                         }
 769                     }
 770                     rmsdav_mat[i][j] = av_tot/weight_tot;
 771                     rmsdav_mat[j][i] = rmsdav_mat[i][j];
 772                     if (rmsdav_mat[i][j] > rmsd_max) rmsd_max=rmsdav_mat[i][j];
 773                 }
 774             }
 775             rmsd_mat=rmsdav_mat;
 776         }
 777
 778         if (bMat) {
 779             fprintf(stderr,"\n%s: Min %f, Max %f, Avg %f\n",
 780                     whatname[ewhat],rmsd_min,rmsd_max,rmsd_avg);
 781             rlo.r = 1; rlo.g = 1; rlo.b = 1;
 782             rhi.r = 0; rhi.g = 0; rhi.b = 0;
 783             if (rmsd_user_max != -1) rmsd_max=rmsd_user_max;
 784             if (rmsd_user_min != -1) rmsd_min=rmsd_user_min;
 785             if ((rmsd_user_max !=  -1) || (rmsd_user_min != -1))
 786                 fprintf(stderr,"Min and Max value set to resp. %f and %f\n",
 787                         rmsd_min,rmsd_max);
 788             sprintf(buf,"%s %s matrix",gn_rms[0],whatname[ewhat]);
 789             write_xpm(opt2FILE("-m",NFILE,fnm,"w"),0,buf,whatlabel[ewhat],
 790                       output_env_get_time_label(oenv),output_env_get_time_label(oenv),tel_mat,tel_mat2,
 791                       axis,axis2, rmsd_mat,rmsd_min,rmsd_max,rlo,rhi,&nlevels);
 792             /* Print the distribution of RMSD values */
 793             if (opt2bSet("-dist",NFILE,fnm))
 794                 low_rmsd_dist(opt2fn("-dist",NFILE,fnm),rmsd_max,tel_mat,rmsd_mat,oenv);
 795
 796             if (bDelta) {
 797                 snew(delta_tot,delta_xsize);
 798                 for(j=0; j<tel_mat-1; j++) {
 799                     for(i=j+1; i<tel_mat; i++) {
 800                         mx=i-j ;
 801                         if (mx < tel_mat/2) {
 802                             if (bDeltaLog)
 803                                 mx=(int)(log(mx)*delta_scalex+0.5);
 804                             my=(int)(rmsd_mat[i][j]*delta_scaley*del_lev+0.5);
 805                             delta_tot[mx] += 1.0;
 806                             if ((rmsd_mat[i][j]>=0) && (rmsd_mat[i][j]<=delta_maxy))
 807                                 delta[mx][my] += 1.0;
 808                         }
 809                     }
 810                 }
 811                 delta_max=0;
 812                 for(i=0; i<delta_xsize; i++) {
 813                     if (delta_tot[i] > 0.0) {
 814                         delta_tot[i]=1.0/delta_tot[i];
 815                         for(j=0; j<=del_lev; j++) {
 816                             delta[i][j] *= delta_tot[i];
 817                             if (delta[i][j] > delta_max)
 818                                 delta_max=delta[i][j];
 819                         }
 820                     }
 821                 }
 822                 fprintf(stderr,"Maximum in delta matrix: %f\n",delta_max);
 823                 snew(del_xaxis,delta_xsize);
 824                 snew(del_yaxis,del_lev+1);
 825                 for (i=0; i<delta_xsize; i++)
 826                     del_xaxis[i]=axis[i]-axis[0];
 827                 for (i=0; i<del_lev+1; i++)
 828                     del_yaxis[i]=delta_maxy*i/del_lev;
 829                 sprintf(buf,"%s %s vs. delta t",gn_rms[0],whatname[ewhat]);
 830                 fp = ffopen("delta.xpm","w");
 831                 write_xpm(fp,0,buf,"density",output_env_get_time_label(oenv),whatlabel[ewhat],
 832                           delta_xsize,del_lev+1,del_xaxis,del_yaxis,
 833                           delta,0.0,delta_max,rlo,rhi,&nlevels);
 834                 ffclose(fp);
 835             }
 836             if (opt2bSet("-bin",NFILE,fnm)) {
 837                 /* NB: File must be binary if we use fwrite */
 838                 fp=ftp2FILE(efDAT,NFILE,fnm,"wb");
 839                 for(i=0;i<tel_mat;i++)
 840                     if(fwrite(rmsd_mat[i],sizeof(**rmsd_mat),tel_mat2,fp) != tel_mat2)
 841                     {
 842                         gmx_fatal(FARGS,"Error writing to output file");
 843                     }
 844                 ffclose(fp);
 845             }
 846         }
 847         if (bBond) {
 848             fprintf(stderr,"\nMin. angle: %f, Max. angle: %f\n",bond_min,bond_max);
 849             if (bond_user_max != -1) bond_max=bond_user_max;
 850             if (bond_user_min != -1) bond_min=bond_user_min;
 851             if ((bond_user_max !=  -1) || (bond_user_min != -1))
 852                 fprintf(stderr,"Bond angle Min and Max set to:\n"
 853                         "Min. angle: %f, Max. angle: %f\n",bond_min,bond_max);
 854             rlo.r = 1; rlo.g = 1; rlo.b = 1;
 855             rhi.r = 0; rhi.g = 0; rhi.b = 0;
 856             sprintf(buf,"%s av. bond angle deviation",gn_rms[0]);
 857             write_xpm(opt2FILE("-bm",NFILE,fnm,"w"),0,buf,"degrees",
 858                       output_env_get_time_label(oenv),output_env_get_time_label(oenv),tel_mat,tel_mat2,
 859                       axis,axis2, bond_mat,bond_min,bond_max,rlo,rhi,&nlevels);
 860         }
 861     }
 862
 863     bAv=opt2bSet("-a",NFILE,fnm);
 864
 865     /* Write the RMSD's to file */
 866     if (!bPrev)
 867         sprintf(buf,"%s",whatxvgname[ewhat]);
 868     else
 869         sprintf(buf,"%s with frame %g %s ago",whatxvgname[ewhat],
 870                 time[prev*freq]-time[0], output_env_get_time_label(oenv));
 871     fp=xvgropen(opt2fn("-o",NFILE,fnm),buf,output_env_get_xvgr_tlabel(oenv),
 872                 whatxvglabel[ewhat], oenv);
 873     if (output_env_get_print_xvgr_codes(oenv))
 874         fprintf(fp,"@ subtitle \"%s%s after %s%s%s\"\n",
 875                 (nrms==1)?"":"of "    , gn_rms[0], fitgraphlabel[efit],
 876                     bFit     ?" to ":""   , bFit?gn_fit:"");
 877     if (nrms != 1)
 878         xvgr_legend(fp,nrms,(const char**)gn_rms,oenv);
 879     for(i=0; (i<teller); i++) {
 880         if ( bSplit && i>0 &&
 881             abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv))<1e-5 )
 882         {
 883             fprintf(fp,"&\n");
 884         }
 885         fprintf(fp,"%12.7f",time[bPrev ? freq*i : i]);
 886         for(j=0; (j<nrms); j++) {
 887             fprintf(fp," %12.7f",rls[j][i]);
 888             if (bAv)
 889                 rlstot+=rls[j][i];
 890         }
 891         fprintf(fp,"\n");
 892     }
 893     ffclose(fp);
 894
 895     if (bMirror) {
 896         /* Write the mirror RMSD's to file */
 897         sprintf(buf,"%s with Mirror",whatxvgname[ewhat]);
 898         sprintf(buf2,"Mirror %s",whatxvglabel[ewhat]);
 899         fp=xvgropen(opt2fn("-mir",NFILE,fnm), buf, output_env_get_xvgr_tlabel(oenv),
 900                     buf2,oenv);
 901         if (nrms == 1) {
 902             if (output_env_get_print_xvgr_codes(oenv))
 903                 fprintf(fp,"@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
 904                         gn_rms[0],gn_fit);
 905         }
 906         else {
 907             if (output_env_get_print_xvgr_codes(oenv))
 908                 fprintf(fp,"@ subtitle \"after lsq fit to mirror %s\"\n",gn_fit);
 909             xvgr_legend(fp,nrms,(const char**)gn_rms,oenv);
 910         }
 911         for(i=0; (i<teller); i++) {
 912             if ( bSplit && i>0 && abs(time[i])<1e-5 )
 913                 fprintf(fp,"&\n");
 914             fprintf(fp,"%12.7f",time[i]);
 915             for(j=0; (j<nrms); j++)
 916                 fprintf(fp," %12.7f",rlsm[j][i]);
 917             fprintf(fp,"\n");
 918         }
 919         ffclose(fp);
 920     }
 921
 922     if (bAv) {
 923         sprintf(buf,"Average %s",whatxvgname[ewhat]);
 924         sprintf(buf2,"Average %s",whatxvglabel[ewhat]);
 925         fp = xvgropen(opt2fn("-a",NFILE,fnm), buf, "Residue", buf2,oenv);
 926         for(j=0; (j<nrms); j++)
 927             fprintf(fp,"%10d  %10g\n",j,rlstot/teller);
 928         ffclose(fp);
 929     }
 930
 931     if (bNorm) {
 932         fp = xvgropen("aver.xvg",gn_rms[0],"Residue",whatxvglabel[ewhat],oenv);
 933         for(j=0; (j<irms[0]); j++)
 934             fprintf(fp,"%10d  %10g\n",j,rlsnorm[j]/teller);
 935         ffclose(fp);
 936     }
 937     do_view(oenv,opt2fn_null("-a",NFILE,fnm),"-graphtype bar");
 938     do_view(oenv,opt2fn("-o",NFILE,fnm),NULL);
 939     do_view(oenv,opt2fn_null("-mir",NFILE,fnm),NULL);
 940     do_view(oenv,opt2fn_null("-m",NFILE,fnm),NULL);
 941     do_view(oenv,opt2fn_null("-bm",NFILE,fnm),NULL);
 942     do_view(oenv,opt2fn_null("-dist",NFILE,fnm),NULL);
 943
 944     thanx(stderr);
 945
 946     return 0;
 947 }