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 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 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 }