src/tools/gmx_cluster.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 #include <math.h>
  39 #include <string.h>
  40 #include <ctype.h>
  41 #include "macros.h"
  42 #include "smalloc.h"
  43 #include "typedefs.h"
  44 #include "copyrite.h"
  45 #include "statutil.h"
  46 #include "tpxio.h"
  47 #include "string2.h"
  48 #include "vec.h"
  49 #include "macros.h"
  50 #include "index.h"
  51 #include "random.h"
  52 #include "pbc.h"
  53 #include "xvgr.h"
  54 #include "futil.h"
  55 #include "matio.h"
  56 #include "eigensolver.h"
  57 #include "cmat.h"
  58 #include "do_fit.h"
  59 #include "trnio.h"
  60 #include "viewit.h"
  61 #include "gmx_ana.h"
  62
  63 /* macro's to print to two file pointers at once (i.e. stderr and log) */
  64 #define lo_ffprintf(fp1,fp2,buf) \
  65    fprintf(fp1,"%s",buf);\
  66    fprintf(fp2,"%s",buf);
  67 /* just print a prepared buffer to fp1 and fp2 */
  68 #define ffprintf(fp1,fp2,buf) { lo_ffprintf(fp1,fp2,buf) }
  69 /* prepare buffer with one argument, then print to fp1 and fp2 */
  70 #define ffprintf1(fp1,fp2,buf,fmt,arg) {\
  71    sprintf(buf,fmt,arg);\
  72    lo_ffprintf(fp1,fp2,buf)\
  73 }
  74 /* prepare buffer with two arguments, then print to fp1 and fp2 */
  75 #define ffprintf2(fp1,fp2,buf,fmt,arg1,arg2) {\
  76    sprintf(buf,fmt,arg1,arg2);\
  77    lo_ffprintf(fp1,fp2,buf)\
  78 }
  79
  80 typedef struct {
  81   int ncl;
  82   int *cl;
  83 } t_clusters;
  84
  85 typedef struct {
  86   int nr;
  87   int *nb;
  88 } t_nnb;
  89
  90 void pr_energy(FILE *fp,real e)
  91 {
  92   fprintf(fp,"Energy: %8.4f\n",e);
  93 }
  94
  95 void cp_index(int nn,int from[],int to[])
  96 {
  97   int i;
  98
  99   for(i=0; (i<nn); i++)
 100     to[i]=from[i];
 101 }
 102
 103 void mc_optimize(FILE *log,t_mat *m,int maxiter,int *seed,real kT)
 104 {
 105   real e[2],ei,ej,efac;
 106   int  *low_index;
 107   int  cur=0;
 108 #define next (1-cur)
 109   int  i,isw,jsw,iisw,jjsw,nn;
 110
 111   fprintf(stderr,"\nDoing Monte Carlo clustering\n");
 112   nn = m->nn;
 113   snew(low_index,nn);
 114   cp_index(nn,m->m_ind,low_index);
 115   if (getenv("TESTMC")) {
 116     e[cur] = mat_energy(m);
 117     pr_energy(log,e[cur]);
 118     fprintf(log,"Doing 1000 random swaps\n");
 119     for(i=0; (i<1000); i++) {
 120       do {
 121         isw = nn*rando(seed);
 122         jsw = nn*rando(seed);
 123       } while ((isw == jsw) || (isw >= nn) || (jsw >= nn));
 124       iisw = m->m_ind[isw];
 125       jjsw = m->m_ind[jsw];
 126       m->m_ind[isw] = jjsw;
 127       m->m_ind[jsw] = iisw;
 128     }
 129   }
 130   e[cur] = mat_energy(m);
 131   pr_energy(log,e[cur]);
 132   for(i=0; (i<maxiter); i++) {
 133     do {
 134       isw = nn*rando(seed);
 135       jsw = nn*rando(seed);
 136     } while ((isw == jsw) || (isw >= nn) || (jsw >= nn));
 137
 138     iisw = m->m_ind[isw];
 139     jjsw = m->m_ind[jsw];
 140     ei   = row_energy(nn,iisw,m->mat[jsw]);
 141     ej   = row_energy(nn,jjsw,m->mat[isw]);
 142
 143     e[next] = e[cur] + (ei+ej-EROW(m,isw)-EROW(m,jsw))/nn;
 144
 145     efac = kT ? exp((e[next]-e[cur])/kT) : -1;
 146     if ((e[next] > e[cur]) || (efac > rando(seed))) {
 147
 148       if (e[next] > e[cur])
 149         cp_index(nn,m->m_ind,low_index);
 150       else
 151         fprintf(log,"Taking uphill step\n");
 152
 153       /* Now swapping rows */
 154       m->m_ind[isw] = jjsw;
 155       m->m_ind[jsw] = iisw;
 156       EROW(m,isw)   = ei;
 157       EROW(m,jsw)   = ej;
 158       cur           = next;
 159       fprintf(log,"Iter: %d Swapped %4d and %4d (now %g)",
 160               i,isw,jsw,mat_energy(m));
 161       pr_energy(log,e[cur]);
 162     }
 163   }
 164   /* Now restore the highest energy index */
 165   cp_index(nn,low_index,m->m_ind);
 166 }
 167
 168 static void calc_dist(int nind,rvec x[],real **d)
 169 {
 170   int     i,j;
 171   real    *xi;
 172   rvec    dx;
 173
 174   for(i=0; (i<nind-1); i++) {
 175     xi=x[i];
 176     for(j=i+1; (j<nind); j++) {
 177       /* Should use pbc_dx when analysing multiple molecueles,
 178        * but the box is not stored for every frame.
 179        */
 180       rvec_sub(xi,x[j],dx);
 181       d[i][j]=norm(dx);
 182     }
 183   }
 184 }
 185
 186 static real rms_dist(int isize,real **d,real **d_r)
 187 {
 188   int  i,j;
 189   real r,r2;
 190
 191   r2=0.0;
 192   for(i=0; (i<isize-1); i++)
 193     for(j=i+1; (j<isize); j++) {
 194       r=d[i][j]-d_r[i][j];
 195       r2+=r*r;
 196     }
 197   r2/=(isize*(isize-1))/2;
 198
 199   return sqrt(r2);
 200 }
 201
 202 static int rms_dist_comp(const void *a,const void *b)
 203 {
 204   t_dist *da,*db;
 205
 206   da = (t_dist *)a;
 207   db = (t_dist *)b;
 208
 209   if (da->dist - db->dist < 0)
 210     return -1;
 211   else if (da->dist - db->dist > 0)
 212     return 1;
 213   return 0;
 214 }
 215
 216 static int clust_id_comp(const void *a,const void *b)
 217 {
 218   t_clustid *da,*db;
 219
 220   da = (t_clustid *)a;
 221   db = (t_clustid *)b;
 222
 223   return da->clust - db->clust;
 224 }
 225
 226 static int nrnb_comp(const void *a, const void *b)
 227 {
 228   t_nnb *da, *db;
 229
 230   da = (t_nnb *)a;
 231   db = (t_nnb *)b;
 232
 233   /* return the b-a, we want highest first */
 234   return db->nr - da->nr;
 235 }
 236
 237 void gather(t_mat *m,real cutoff,t_clusters *clust)
 238 {
 239   t_clustid *c;
 240   t_dist    *d;
 241   int       i,j,k,nn,cid,n1,diff;
 242   gmx_bool      bChange;
 243
 244   /* First we sort the entries in the RMSD matrix */
 245   n1 = m->nn;
 246   nn = ((n1-1)*n1)/2;
 247   snew(d,nn);
 248   for(i=k=0; (i<n1); i++)
 249     for(j=i+1; (j<n1); j++,k++) {
 250       d[k].i    = i;
 251       d[k].j    = j;
 252       d[k].dist = m->mat[i][j];
 253     }
 254   if (k != nn)
 255     gmx_incons("gather algortihm");
 256   qsort(d,nn,sizeof(d[0]),rms_dist_comp);
 257
 258   /* Now we make a cluster index for all of the conformations */
 259   c = new_clustid(n1);
 260
 261   /* Now we check the closest structures, and equalize their cluster numbers */
 262   fprintf(stderr,"Linking structures ");
 263   do {
 264     fprintf(stderr,"*");
 265     bChange=FALSE;
 266     for(k=0; (k<nn) && (d[k].dist < cutoff); k++) {
 267       diff = c[d[k].j].clust - c[d[k].i].clust;
 268       if (diff) {
 269         bChange = TRUE;
 270         if (diff > 0)
 271           c[d[k].j].clust = c[d[k].i].clust;
 272         else
 273           c[d[k].i].clust = c[d[k].j].clust;
 274       }
 275     }
 276   } while (bChange);
 277   fprintf(stderr,"\nSorting and renumbering clusters\n");
 278   /* Sort on cluster number */
 279   qsort(c,n1,sizeof(c[0]),clust_id_comp);
 280
 281   /* Renumber clusters */
 282   cid = 1;
 283   for(k=1; k<n1; k++) {
 284     if (c[k].clust != c[k-1].clust) {
 285       c[k-1].clust = cid;
 286       cid ++;
 287     } else
 288       c[k-1].clust = cid;
 289   }
 290   c[k-1].clust = cid;
 291   if (debug)
 292     for(k=0; (k<n1); k++)
 293       fprintf(debug,"Cluster index for conformation %d: %d\n",
 294               c[k].conf,c[k].clust);
 295   clust->ncl = cid;
 296   for(k=0; k<n1; k++)
 297     clust->cl[c[k].conf] = c[k].clust;
 298
 299   sfree(c);
 300   sfree(d);
 301 }
 302
 303 gmx_bool jp_same(int **nnb,int i,int j,int P)
 304 {
 305   gmx_bool bIn;
 306   int  k,ii,jj,pp;
 307
 308   bIn = FALSE;
 309   for(k=0; nnb[i][k]>=0; k++)
 310     bIn = bIn || (nnb[i][k] == j);
 311   if (!bIn)
 312     return FALSE;
 313
 314   bIn = FALSE;
 315   for(k=0; nnb[j][k]>=0; k++)
 316     bIn = bIn || (nnb[j][k] == i);
 317   if (!bIn)
 318     return FALSE;
 319
 320   pp=0;
 321   for(ii=0; nnb[i][ii]>=0; ii++)
 322     for(jj=0; nnb[j][jj]>=0; jj++)
 323       if ((nnb[i][ii] == nnb[j][jj]) && (nnb[i][ii] != -1))
 324         pp++;
 325
 326   return (pp >= P);
 327 }
 328
 329 static void jarvis_patrick(int n1,real **mat,int M,int P,
 330                            real rmsdcut,t_clusters *clust)
 331 {
 332   t_dist    *row;
 333   t_clustid *c;
 334   int       **nnb;
 335   int       i,j,k,cid,diff,max;
 336   gmx_bool      bChange;
 337   real      **mcpy=NULL;
 338
 339   if (rmsdcut < 0)
 340     rmsdcut = 10000;
 341
 342   /* First we sort the entries in the RMSD matrix row by row.
 343    * This gives us the nearest neighbor list.
 344    */
 345   snew(nnb,n1);
 346   snew(row,n1);
 347   for(i=0; (i<n1); i++) {
 348     for(j=0; (j<n1); j++) {
 349       row[j].j    = j;
 350       row[j].dist = mat[i][j];
 351     }
 352     qsort(row,n1,sizeof(row[0]),rms_dist_comp);
 353     if (M>0) {
 354       /* Put the M nearest neighbors in the list */
 355       snew(nnb[i],M+1);
 356       for(j=k=0; (k<M) && (j<n1) && (mat[i][row[j].j] < rmsdcut); j++)
 357         if (row[j].j  != i) {
 358           nnb[i][k]  = row[j].j;
 359           k++;
 360         }
 361       nnb[i][k] = -1;
 362     } else {
 363       /* Put all neighbors nearer than rmsdcut in the list */
 364       max=0;
 365       k=0;
 366       for(j=0; (j<n1) && (mat[i][row[j].j] < rmsdcut); j++)
 367         if (row[j].j != i) {
 368           if (k >= max) {
 369             max += 10;
 370             srenew(nnb[i],max);
 371           }
 372           nnb[i][k] = row[j].j;
 373           k++;
 374         }
 375       if (k == max)
 376         srenew(nnb[i],max+1);
 377       nnb[i][k] = -1;
 378     }
 379   }
 380   sfree(row);
 381   if (debug) {
 382     fprintf(debug,"Nearest neighborlist. M = %d, P = %d\n",M,P);
 383     for(i=0; (i<n1); i++) {
 384       fprintf(debug,"i:%5d nbs:",i);
 385       for(j=0; nnb[i][j]>=0; j++)
 386         fprintf(debug,"%5d[%5.3f]",nnb[i][j],mat[i][nnb[i][j]]);
 387       fprintf(debug,"\n");
 388     }
 389   }
 390
 391   c = new_clustid(n1);
 392   fprintf(stderr,"Linking structures ");
 393   /* Use mcpy for temporary storage of gmx_booleans */
 394   mcpy = mk_matrix(n1,n1,FALSE);
 395   for(i=0; i<n1; i++)
 396     for(j=i+1; j<n1; j++)
 397       mcpy[i][j] = jp_same(nnb,i,j,P);
 398   do {
 399     fprintf(stderr,"*");
 400     bChange=FALSE;
 401     for(i=0; i<n1; i++) {
 402       for(j=i+1; j<n1; j++)
 403         if (mcpy[i][j]) {
 404           diff = c[j].clust - c[i].clust;
 405           if (diff) {
 406             bChange = TRUE;
 407             if (diff > 0)
 408               c[j].clust = c[i].clust;
 409             else
 410               c[i].clust = c[j].clust;
 411           }
 412         }
 413     }
 414   } while (bChange);
 415
 416   fprintf(stderr,"\nSorting and renumbering clusters\n");
 417   /* Sort on cluster number */
 418   qsort(c,n1,sizeof(c[0]),clust_id_comp);
 419
 420   /* Renumber clusters */
 421   cid = 1;
 422   for(k=1; k<n1; k++) {
 423     if (c[k].clust != c[k-1].clust) {
 424       c[k-1].clust = cid;
 425       cid ++;
 426     } else
 427       c[k-1].clust = cid;
 428   }
 429   c[k-1].clust = cid;
 430   clust->ncl = cid;
 431   for(k=0; k<n1; k++)
 432     clust->cl[c[k].conf] = c[k].clust;
 433   if (debug)
 434     for(k=0; (k<n1); k++)
 435       fprintf(debug,"Cluster index for conformation %d: %d\n",
 436               c[k].conf,c[k].clust);
 437
 438 /* Again, I don't see the point in this... (AF) */
 439 /*   for(i=0; (i<n1); i++) { */
 440 /*     for(j=0; (j<n1); j++) */
 441 /*       mcpy[c[i].conf][c[j].conf] = mat[i][j]; */
 442 /*   } */
 443 /*   for(i=0; (i<n1); i++) { */
 444 /*     for(j=0; (j<n1); j++) */
 445 /*       mat[i][j] = mcpy[i][j]; */
 446 /*   } */
 447   done_matrix(n1,&mcpy);
 448
 449   sfree(c);
 450   for(i=0; (i<n1); i++)
 451     sfree(nnb[i]);
 452   sfree(nnb);
 453 }
 454
 455 static void dump_nnb (FILE *fp, const char *title, int n1, t_nnb *nnb)
 456 {
 457   int i,j;
 458
 459   /* dump neighbor list */
 460   fprintf(fp,"%s",title);
 461   for(i=0; (i<n1); i++) {
 462     fprintf(fp,"i:%5d #:%5d nbs:",i,nnb[i].nr);
 463     for(j=0; j<nnb[i].nr; j++)
 464       fprintf(fp,"%5d",nnb[i].nb[j]);
 465     fprintf(fp,"\n");
 466   }
 467 }
 468
 469 static void gromos(int n1, real **mat, real rmsdcut, t_clusters *clust)
 470 {
 471   t_dist *row;
 472   t_nnb  *nnb;
 473   int    i,j,k,j1,max;
 474
 475   /* Put all neighbors nearer than rmsdcut in the list */
 476   fprintf(stderr,"Making list of neighbors within cutoff ");
 477   snew(nnb,n1);
 478   snew(row,n1);
 479   for(i=0; (i<n1); i++) {
 480     max=0;
 481     k=0;
 482     /* put all neighbors within cut-off in list */
 483     for(j=0; j<n1; j++)
 484       if (mat[i][j] < rmsdcut) {
 485         if (k >= max) {
 486           max += 10;
 487           srenew(nnb[i].nb,max);
 488         }
 489         nnb[i].nb[k] = j;
 490         k++;
 491       }
 492     /* store nr of neighbors, we'll need that */
 493     nnb[i].nr = k;
 494     if (i%(1+n1/100)==0) fprintf(stderr,"%3d%%\b\b\b\b",(i*100+1)/n1);
 495   }
 496   fprintf(stderr,"%3d%%\n",100);
 497   sfree(row);
 498
 499   /* sort neighbor list on number of neighbors, largest first */
 500   qsort(nnb,n1,sizeof(nnb[0]),nrnb_comp);
 501
 502   if (debug) dump_nnb(debug, "Nearest neighborlist after sort.\n", n1, nnb);
 503
 504   /* turn first structure with all its neighbors (largest) into cluster
 505      remove them from pool of structures and repeat for all remaining */
 506   fprintf(stderr,"Finding clusters %4d", 0);
 507   /* cluster id's start at 1: */
 508   k=1;
 509   while(nnb[0].nr) {
 510     /* set cluster id (k) for first item in neighborlist */
 511     for (j=0; j<nnb[0].nr; j++)
 512       clust->cl[nnb[0].nb[j]] = k;
 513     /* mark as done */
 514     nnb[0].nr=0;
 515     sfree(nnb[0].nb);
 516
 517     /* adjust number of neighbors for others, taking removals into account: */
 518     for(i=1; i<n1 && nnb[i].nr; i++) {
 519       j1=0;
 520       for(j=0; j<nnb[i].nr; j++)
 521         /* if this neighbor wasn't removed */
 522         if ( clust->cl[nnb[i].nb[j]] == 0 ) {
 523           /* shift the rest (j1<=j) */
 524           nnb[i].nb[j1]=nnb[i].nb[j];
 525           /* next */
 526           j1++;
 527         }
 528       /* now j1 is the new number of neighbors */
 529       nnb[i].nr=j1;
 530     }
 531     /* sort again on nnb[].nr, because we have new # neighbors: */
 532     /* but we only need to sort upto i, i.e. when nnb[].nr>0 */
 533     qsort(nnb,i,sizeof(nnb[0]),nrnb_comp);
 534
 535     fprintf(stderr,"\b\b\b\b%4d",k);
 536     /* new cluster id */
 537     k++;
 538   }
 539   fprintf(stderr,"\n");
 540   sfree(nnb);
 541   if (debug) {
 542     fprintf(debug,"Clusters (%d):\n", k);
 543     for(i=0; i<n1; i++)
 544       fprintf(debug," %3d", clust->cl[i]);
 545     fprintf(debug,"\n");
 546   }
 547
 548   clust->ncl=k-1;
 549 }
 550
 551 rvec **read_whole_trj(const char *fn,int isize,atom_id index[],int skip,
 552                       int *nframe, real **time,const output_env_t oenv)
 553 {
 554   rvec   **xx,*x;
 555   matrix box;
 556   real   t;
 557   int    i,i0,j,max_nf;
 558   int    natom;
 559   t_trxstatus *status;
 560
 561
 562   max_nf = 0;
 563   xx     = NULL;
 564   *time  = NULL;
 565   natom = read_first_x(oenv,&status,fn,&t,&x,box);
 566   i  = 0;
 567   i0 = 0;
 568   do {
 569     if (i0 >= max_nf) {
 570       max_nf += 10;
 571       srenew(xx,max_nf);
 572       srenew(*time,max_nf);
 573     }
 574     if ((i % skip) == 0) {
 575       snew(xx[i0],isize);
 576       /* Store only the interesting atoms */
 577       for(j=0; (j<isize); j++)
 578         copy_rvec(x[index[j]],xx[i0][j]);
 579       (*time)[i0] = t;
 580       i0 ++;
 581     }
 582     i++;
 583   } while (read_next_x(oenv,status,&t,natom,x,box));
 584   fprintf(stderr,"Allocated %lu bytes for frames\n",
 585           (unsigned long) (max_nf*isize*sizeof(**xx)));
 586   fprintf(stderr,"Read %d frames from trajectory %s\n",i0,fn);
 587   *nframe = i0;
 588   sfree(x);
 589
 590   return xx;
 591 }
 592
 593 static int plot_clusters(int nf, real **mat, t_clusters *clust,
 594                          int nlevels, int minstruct)
 595 {
 596   int i,j,ncluster,ci;
 597   int *cl_id,*nstruct,*strind;
 598
 599   snew(cl_id,nf);
 600   snew(nstruct,nf);
 601   snew(strind,nf);
 602   for(i=0; i<nf; i++) {
 603     strind[i] = 0;
 604     cl_id[i]  = clust->cl[i];
 605     nstruct[cl_id[i]]++;
 606   }
 607   ncluster = 0;
 608   for(i=0; i<nf; i++) {
 609     if (nstruct[i] >= minstruct) {
 610       ncluster++;
 611       for(j=0; (j<nf); j++)
 612         if (cl_id[j] == i)
 613           strind[j] = ncluster;
 614     }
 615   }
 616   ncluster++;
 617   fprintf(stderr,"There are %d clusters with at least %d conformations\n",
 618             ncluster,minstruct);
 619
 620   for(i=0; (i<nf); i++) {
 621     ci = cl_id[i];
 622     for(j=0; j<i; j++)
 623       if ((ci == cl_id[j]) && (nstruct[ci] >= minstruct)) {
 624         /* color different clusters with different colors, as long as
 625            we don't run out of colors */
 626         mat[i][j] = strind[i];
 627       }
 628       else
 629         mat[i][j] = 0;
 630   }
 631   sfree(strind);
 632   sfree(nstruct);
 633   sfree(cl_id);
 634
 635   return ncluster;
 636 }
 637
 638 static void mark_clusters(int nf, real **mat, real val, t_clusters *clust)
 639 {
 640   int i,j,v;
 641
 642   for(i=0; i<nf; i++)
 643     for(j=0; j<i; j++)
 644       if (clust->cl[i] == clust->cl[j])
 645         mat[i][j] = val;
 646       else
 647         mat[i][j] = 0;
 648 }
 649
 650 static char *parse_filename(const char *fn, int maxnr)
 651 {
 652   int i;
 653   char *fnout;
 654   const char *ext;
 655   char buf[STRLEN];
 656
 657   if (strchr(fn,'%'))
 658     gmx_fatal(FARGS,"will not number filename %s containing '%c'",fn,'%');
 659   /* number of digits needed in numbering */
 660   i = (int)(log(maxnr)/log(10)) + 1;
 661   /* split fn and ext */
 662   ext = strrchr(fn, '.');
 663   if (!ext)
 664     gmx_fatal(FARGS,"cannot separate extension in filename %s",fn);
 665   ext++;
 666   /* insert e.g. '%03d' between fn and ext */
 667   sprintf(buf,"%s%%0%dd.%s",fn,i,ext);
 668   snew(fnout,strlen(buf)+1);
 669   strcpy(fnout, buf);
 670
 671   return fnout;
 672 }
 673
 674 static void ana_trans(t_clusters *clust, int nf,
 675                       const char *transfn, const char *ntransfn, FILE *log,
 676                       t_rgb rlo,t_rgb rhi,const output_env_t oenv)
 677 {
 678   FILE *fp;
 679   real **trans,*axis;
 680   int  *ntrans;
 681   int  i,ntranst,maxtrans;
 682   char buf[STRLEN];
 683
 684   snew(ntrans,clust->ncl);
 685   snew(trans,clust->ncl);
 686   snew(axis,clust->ncl);
 687   for(i=0; i<clust->ncl; i++) {
 688     axis[i]=i+1;
 689     snew(trans[i],clust->ncl);
 690   }
 691   ntranst=0;
 692   maxtrans=0;
 693   for(i=1; i<nf; i++)
 694     if(clust->cl[i] != clust->cl[i-1]) {
 695       ntranst++;
 696       ntrans[clust->cl[i-1]-1]++;
 697       ntrans[clust->cl[i]-1]++;
 698       trans[clust->cl[i-1]-1][clust->cl[i]-1]++;
 699       maxtrans = max(maxtrans, trans[clust->cl[i]-1][clust->cl[i-1]-1]);
 700     }
 701   ffprintf2(stderr,log,buf,"Counted %d transitions in total, "
 702             "max %d between two specific clusters\n",ntranst,maxtrans);
 703   if (transfn) {
 704     fp=ffopen(transfn,"w");
 705     i = min(maxtrans+1, 80);
 706     write_xpm(fp,0,"Cluster Transitions","# transitions",
 707               "from cluster","to cluster",
 708               clust->ncl, clust->ncl, axis, axis, trans,
 709               0, maxtrans, rlo, rhi, &i);
 710     ffclose(fp);
 711   }
 712   if (ntransfn) {
 713     fp=xvgropen(ntransfn,"Cluster Transitions","Cluster #","# transitions",
 714                 oenv);
 715     for(i=0; i<clust->ncl; i++)
 716       fprintf(fp,"%5d %5d\n",i+1,ntrans[i]);
 717     ffclose(fp);
 718   }
 719   sfree(ntrans);
 720   for(i=0; i<clust->ncl; i++)
 721     sfree(trans[i]);
 722   sfree(trans);
 723   sfree(axis);
 724 }
 725
 726 static void analyze_clusters(int nf, t_clusters *clust, real **rmsd,
 727                              int natom, t_atoms *atoms, rvec *xtps,
 728                              real *mass, rvec **xx, real *time,
 729                              int ifsize, atom_id *fitidx,
 730                              int iosize, atom_id *outidx,
 731                              const char *trxfn, const char *sizefn,
 732                              const char *transfn, const char *ntransfn,
 733                              const char *clustidfn, gmx_bool bAverage,
 734                              int write_ncl, int write_nst, real rmsmin,
 735                              gmx_bool bFit, FILE *log,t_rgb rlo,t_rgb rhi,
 736                              const output_env_t oenv)
 737 {
 738   FILE *fp=NULL;
 739   char buf[STRLEN],buf1[40],buf2[40],buf3[40],*trxsfn;
 740   t_trxstatus *trxout=NULL;
 741   t_trxstatus *trxsout=NULL;
 742   int  i,i1,cl,nstr,*structure,first=0,midstr;
 743   gmx_bool *bWrite=NULL;
 744   real r,clrmsd,midrmsd;
 745   rvec *xav=NULL;
 746   matrix zerobox;
 747
 748   clear_mat(zerobox);
 749
 750   ffprintf1(stderr,log,buf,"\nFound %d clusters\n\n",clust->ncl);
 751   trxsfn=NULL;
 752   if (trxfn) {
 753     /* do we write all structures? */
 754     if (write_ncl) {
 755       trxsfn = parse_filename(trxfn, max(write_ncl,clust->ncl));
 756       snew(bWrite,nf);
 757     }
 758     ffprintf2(stderr,log,buf,"Writing %s structure for each cluster to %s\n",
 759               bAverage ? "average" : "middle", trxfn);
 760     if (write_ncl) {
 761       /* find out what we want to tell the user:
 762          Writing [all structures|structures with rmsd > %g] for
 763          {all|first %d} clusters {with more than %d structures} to %s     */
 764       if (rmsmin>0.0)
 765         sprintf(buf1,"structures with rmsd > %g",rmsmin);
 766       else
 767         sprintf(buf1,"all structures");
 768       buf2[0]=buf3[0]='\0';
 769       if (write_ncl>=clust->ncl) {
 770         if (write_nst==0)
 771           sprintf(buf2,"all ");
 772       } else
 773         sprintf(buf2,"the first %d ",write_ncl);
 774       if (write_nst)
 775         sprintf(buf3," with more than %d structures",write_nst);
 776       sprintf(buf,"Writing %s for %sclusters%s to %s\n",buf1,buf2,buf3,trxsfn);
 777       ffprintf(stderr,log,buf);
 778     }
 779
 780     /* Prepare a reference structure for the orientation of the clusters  */
 781     if (bFit)
 782         reset_x(ifsize,fitidx,natom,NULL,xtps,mass);
 783     trxout = open_trx(trxfn,"w");
 784     /* Calculate the average structure in each cluster,               *
 785      * all structures are fitted to the first struture of the cluster */
 786     snew(xav,natom);
 787   }
 788
 789   if (transfn || ntransfn)
 790     ana_trans(clust, nf, transfn, ntransfn, log,rlo,rhi,oenv);
 791
 792   if (clustidfn) {
 793     fp=xvgropen(clustidfn,"Clusters",output_env_get_xvgr_tlabel(oenv),"Cluster #",oenv);
 794     fprintf(fp,"@    s0 symbol 2\n");
 795     fprintf(fp,"@    s0 symbol size 0.2\n");
 796     fprintf(fp,"@    s0 linestyle 0\n");
 797     for(i=0; i<nf; i++)
 798       fprintf(fp,"%8g %8d\n",time[i],clust->cl[i]);
 799     ffclose(fp);
 800   }
 801   if (sizefn) {
 802     fp=xvgropen(sizefn,"Cluster Sizes","Cluster #","# Structures",oenv);
 803     fprintf(fp,"@g%d type %s\n",0,"bar");
 804   }
 805   snew(structure,nf);
 806   fprintf(log,"\n%3s | %3s  %4s | %6s %4s | cluster members\n",
 807           "cl.","#st","rmsd","middle","rmsd");
 808   for(cl=1; cl<=clust->ncl; cl++) {
 809     /* prepare structures (fit, middle, average) */
 810     if (xav)
 811       for(i=0; i<natom;i++)
 812         clear_rvec(xav[i]);
 813     nstr=0;
 814     for(i1=0; i1<nf; i1++)
 815       if (clust->cl[i1] == cl) {
 816         structure[nstr] = i1;
 817         nstr++;
 818         if (trxfn && (bAverage || write_ncl) ) {
 819           if (bFit)
 820           reset_x(ifsize,fitidx,natom,NULL,xx[i1],mass);
 821           if (nstr == 1)
 822             first = i1;
 823           else if (bFit)
 824             do_fit(natom,mass,xx[first],xx[i1]);
 825           if (xav)
 826             for(i=0; i<natom; i++)
 827               rvec_inc(xav[i],xx[i1][i]);
 828         }
 829       }
 830     if (sizefn)
 831       fprintf(fp,"%8d %8d\n",cl,nstr);
 832     clrmsd = 0;
 833     midstr = 0;
 834     midrmsd = 10000;
 835     for(i1=0; i1<nstr; i1++) {
 836       r = 0;
 837       if (nstr > 1) {
 838         for(i=0; i<nstr; i++)
 839           if (i < i1)
 840             r += rmsd[structure[i]][structure[i1]];
 841           else
 842             r += rmsd[structure[i1]][structure[i]];
 843         r /= (nstr - 1);
 844       }
 845       if ( r < midrmsd ) {
 846         midstr = structure[i1];
 847         midrmsd = r;
 848       }
 849       clrmsd += r;
 850     }
 851     clrmsd /= nstr;
 852
 853     /* dump cluster info to logfile */
 854     if (nstr > 1) {
 855       sprintf(buf1,"%6.3f",clrmsd);
 856       if (buf1[0] == '0')
 857         buf1[0] = ' ';
 858       sprintf(buf2,"%5.3f",midrmsd);
 859       if (buf2[0] == '0')
 860         buf2[0] = ' ';
 861     } else {
 862       sprintf(buf1,"%5s","");
 863       sprintf(buf2,"%5s","");
 864     }
 865     fprintf(log,"%3d | %3d %s | %6g%s |",cl,nstr,buf1,time[midstr],buf2);
 866     for(i=0; i<nstr; i++) {
 867       if ((i % 7 == 0) && i)
 868         sprintf(buf,"\n%3s | %3s  %4s | %6s %4s |","","","","","");
 869       else
 870         buf[0] = '\0';
 871       i1 = structure[i];
 872       fprintf(log,"%s %6g",buf,time[i1]);
 873     }
 874     fprintf(log,"\n");
 875
 876     /* write structures to trajectory file(s) */
 877     if (trxfn) {
 878       if (write_ncl)
 879         for(i=0; i<nstr; i++)
 880           bWrite[i]=FALSE;
 881       if ( cl < write_ncl+1 && nstr > write_nst ) {
 882         /* Dump all structures for this cluster */
 883         /* generate numbered filename (there is a %d in trxfn!) */
 884         sprintf(buf,trxsfn,cl);
 885         trxsout = open_trx(buf,"w");
 886         for(i=0; i<nstr; i++) {
 887           bWrite[i] = TRUE;
 888           if (rmsmin>0.0)
 889             for(i1=0; i1<i && bWrite[i]; i1++)
 890               if (bWrite[i1])
 891                 bWrite[i] = rmsd[structure[i1]][structure[i]] > rmsmin;
 892           if (bWrite[i])
 893             write_trx(trxsout,iosize,outidx,atoms,i,time[structure[i]],zerobox,
 894                       xx[structure[i]],NULL,NULL);
 895         }
 896         close_trx(trxsout);
 897       }
 898       /* Dump the average structure for this cluster */
 899       if (bAverage) {
 900         for(i=0; i<natom; i++)
 901           svmul(1.0/nstr,xav[i],xav[i]);
 902       } else {
 903         for(i=0; i<natom; i++)
 904           copy_rvec(xx[midstr][i],xav[i]);
 905         if (bFit)
 906         reset_x(ifsize,fitidx,natom,NULL,xav,mass);
 907       }
 908       if (bFit)
 909         do_fit(natom,mass,xtps,xav);
 910       r = cl;
 911       write_trx(trxout,iosize,outidx,atoms,cl,time[midstr],zerobox,xav,NULL,NULL);
 912     }
 913   }
 914   /* clean up */
 915   if (trxfn) {
 916     close_trx(trxout);
 917     sfree(xav);
 918     if (write_ncl)
 919       sfree(bWrite);
 920   }
 921   sfree(structure);
 922   if (trxsfn)
 923     sfree(trxsfn);
 924 }
 925
 926 static void convert_mat(t_matrix *mat,t_mat *rms)
 927 {
 928   int i,j;
 929
 930   rms->n1 = mat->nx;
 931   matrix2real(mat,rms->mat);
 932   /* free input xpm matrix data */
 933   for(i=0; i<mat->nx; i++)
 934     sfree(mat->matrix[i]);
 935   sfree(mat->matrix);
 936
 937   for(i=0; i<mat->nx; i++)
 938     for(j=i; j<mat->nx; j++) {
 939       rms->sumrms += rms->mat[i][j];
 940       rms->maxrms = max(rms->maxrms, rms->mat[i][j]);
 941       if (j!=i)
 942         rms->minrms = min(rms->minrms, rms->mat[i][j]);
 943     }
 944   rms->nn = mat->nx;
 945 }
 946
 947 int gmx_cluster(int argc,char *argv[])
 948 {
 949   const char *desc[] = {
 950     "g_cluster can cluster structures with several different methods.",
 951     "Distances between structures can be determined from a trajectory",
 952     "or read from an XPM matrix file with the [TT]-dm[tt] option.",
 953     "RMS deviation after fitting or RMS deviation of atom-pair distances",
 954     "can be used to define the distance between structures.[PAR]",
 955
 956     "single linkage: add a structure to a cluster when its distance to any",
 957     "element of the cluster is less than [TT]cutoff[tt].[PAR]",
 958
 959     "Jarvis Patrick: add a structure to a cluster when this structure",
 960     "and a structure in the cluster have each other as neighbors and",
 961     "they have a least [TT]P[tt] neighbors in common. The neighbors",
 962     "of a structure are the M closest structures or all structures within",
 963     "[TT]cutoff[tt].[PAR]",
 964
 965     "Monte Carlo: reorder the RMSD matrix using Monte Carlo.[PAR]",
 966
 967     "diagonalization: diagonalize the RMSD matrix.[PAR]",
 968
 969     "gromos: use algorithm as described in Daura [IT]et al.[it]",
 970     "([IT]Angew. Chem. Int. Ed.[it] [BB]1999[bb], [IT]38[it], pp 236-240).",
 971     "Count number of neighbors using cut-off, take structure with",
 972     "largest number of neighbors with all its neighbors as cluster",
 973     "and eleminate it from the pool of clusters. Repeat for remaining",
 974     "structures in pool.[PAR]",
 975
 976     "When the clustering algorithm assigns each structure to exactly one",
 977     "cluster (single linkage, Jarvis Patrick and gromos) and a trajectory",
 978     "file is supplied, the structure with",
 979     "the smallest average distance to the others or the average structure",
 980     "or all structures for each cluster will be written to a trajectory",
 981     "file. When writing all structures, separate numbered files are made",
 982     "for each cluster.[PAR]",
 983
 984     "Two output files are always written:[BR]",
 985     "[TT]-o[tt] writes the RMSD values in the upper left half of the matrix",
 986     "and a graphical depiction of the clusters in the lower right half",
 987     "When [TT]-minstruct[tt] = 1 the graphical depiction is black",
 988     "when two structures are in the same cluster.",
 989     "When [TT]-minstruct[tt] > 1 different colors will be used for each",
 990     "cluster.[BR]",
 991     "[TT]-g[tt] writes information on the options used and a detailed list",
 992     "of all clusters and their members.[PAR]",
 993
 994     "Additionally, a number of optional output files can be written:[BR]",
 995     "[TT]-dist[tt] writes the RMSD distribution.[BR]",
 996     "[TT]-ev[tt] writes the eigenvectors of the RMSD matrix",
 997     "diagonalization.[BR]",
 998     "[TT]-sz[tt] writes the cluster sizes.[BR]",
 999     "[TT]-tr[tt] writes a matrix of the number transitions between",
1000     "cluster pairs.[BR]",
1001     "[TT]-ntr[tt] writes the total number of transitions to or from",
1002     "each cluster.[BR]",
1003     "[TT]-clid[tt] writes the cluster number as a function of time.[BR]",
1004     "[TT]-cl[tt] writes average (with option [TT]-av[tt]) or central",
1005     "structure of each cluster or writes numbered files with cluster members",
1006     "for a selected set of clusters (with option [TT]-wcl[tt], depends on",
1007     "[TT]-nst[tt] and [TT]-rmsmin[tt]).[BR]",
1008   };
1009
1010   FILE         *fp,*log;
1011   int          i,i1,i2,j,nf,nrms;
1012
1013   matrix       box;
1014   rvec         *xtps,*usextps,*x1,**xx=NULL;
1015   const char   *fn,*trx_out_fn;
1016   t_clusters   clust;
1017   t_mat        *rms;
1018   real         *eigval;
1019   t_topology   top;
1020   int          ePBC;
1021   t_atoms      useatoms;
1022   t_matrix     *readmat=NULL;
1023   real         *tmp;
1024
1025   int      isize=0,ifsize=0,iosize=0;
1026   atom_id  *index=NULL, *fitidx, *outidx;
1027   char     *grpname;
1028   real     rmsd,**d1,**d2,*time=NULL,time_invfac,*mass=NULL;
1029   char     buf[STRLEN],buf1[80],title[STRLEN];
1030   gmx_bool     bAnalyze,bUseRmsdCut,bJP_RMSD=FALSE,bReadMat,bReadTraj;
1031
1032   int method,ncluster=0;
1033   static const char *methodname[] = {
1034     NULL, "linkage", "jarvis-patrick","monte-carlo",
1035     "diagonalization", "gromos", NULL
1036   };
1037   enum { m_null, m_linkage, m_jarvis_patrick,
1038          m_monte_carlo, m_diagonalize, m_gromos, m_nr };
1039   /* Set colors for plotting: white = zero RMS, black = maximum */
1040   static t_rgb rlo_top = { 1.0, 1.0, 1.0 };
1041   static t_rgb rhi_top = { 0.0, 0.0, 0.0 };
1042   static t_rgb rlo_bot = { 1.0, 1.0, 1.0 };
1043   static t_rgb rhi_bot = { 0.0, 0.0, 1.0 };
1044   static int  nlevels=40,skip=1;
1045   static real scalemax=-1.0,rmsdcut=0.1,rmsmin=0.0;
1046   static gmx_bool bRMSdist=FALSE,bBinary=FALSE,bAverage=FALSE,bFit=TRUE;
1047   static int  niter=10000,seed=1993,write_ncl=0,write_nst=1,minstruct=1;
1048   static real kT=1e-3;
1049   static int  M=10,P=3;
1050   output_env_t oenv;
1051   t_pargs pa[] = {
1052     { "-dista", FALSE, etBOOL, {&bRMSdist},
1053       "Use RMSD of distances instead of RMS deviation" },
1054     { "-nlevels",FALSE,etINT,  {&nlevels},
1055       "Discretize RMSD matrix in # levels" },
1056     { "-cutoff",FALSE, etREAL, {&rmsdcut},
1057       "RMSD cut-off (nm) for two structures to be neighbor" },
1058     { "-fit",   FALSE, etBOOL, {&bFit},
1059       "Use least squares fitting before RMSD calculation" },
1060     { "-max",   FALSE, etREAL, {&scalemax},
1061       "Maximum level in RMSD matrix" },
1062     { "-skip",  FALSE, etINT,  {&skip},
1063       "Only analyze every nr-th frame" },
1064     { "-av",    FALSE, etBOOL, {&bAverage},
1065       "Write average iso middle structure for each cluster" },
1066     { "-wcl",   FALSE, etINT,  {&write_ncl},
1067       "Write all structures for first # clusters to numbered files" },
1068     { "-nst",   FALSE, etINT,  {&write_nst},
1069       "Only write all structures if more than # per cluster" },
1070     { "-rmsmin",FALSE, etREAL, {&rmsmin},
1071       "minimum rms difference with rest of cluster for writing structures" },
1072     { "-method",FALSE, etENUM, {methodname},
1073       "Method for cluster determination" },
1074     { "-minstruct", FALSE, etINT, {&minstruct},
1075       "Minimum number of structures in cluster for coloring in the xpm file" },
1076     { "-binary",FALSE, etBOOL, {&bBinary},
1077       "Treat the RMSD matrix as consisting of 0 and 1, where the cut-off "
1078       "is given by -cutoff" },
1079     { "-M",     FALSE, etINT,  {&M},
1080       "Number of nearest neighbors considered for Jarvis-Patrick algorithm, "
1081       "0 is use cutoff" },
1082     { "-P",     FALSE, etINT,  {&P},
1083       "Number of identical nearest neighbors required to form a cluster" },
1084     { "-seed",  FALSE, etINT,  {&seed},
1085       "Random number seed for Monte Carlo clustering algorithm" },
1086     { "-niter", FALSE, etINT,  {&niter},
1087       "Number of iterations for MC" },
1088     { "-kT",    FALSE, etREAL, {&kT},
1089       "Boltzmann weighting factor for Monte Carlo optimization "
1090       "(zero turns off uphill steps)" }
1091   };
1092   t_filenm fnm[] = {
1093     { efTRX, "-f",     NULL,        ffOPTRD },
1094     { efTPS, "-s",     NULL,        ffOPTRD },
1095     { efNDX, NULL,     NULL,        ffOPTRD },
1096     { efXPM, "-dm",   "rmsd",       ffOPTRD },
1097     { efXPM, "-o",    "rmsd-clust", ffWRITE },
1098     { efLOG, "-g",    "cluster",    ffWRITE },
1099     { efXVG, "-dist", "rmsd-dist",  ffOPTWR },
1100     { efXVG, "-ev",   "rmsd-eig",   ffOPTWR },
1101     { efXVG, "-sz",   "clust-size", ffOPTWR},
1102     { efXPM, "-tr",   "clust-trans",ffOPTWR},
1103     { efXVG, "-ntr",  "clust-trans",ffOPTWR},
1104     { efXVG, "-clid", "clust-id.xvg",ffOPTWR},
1105     { efTRX, "-cl",   "clusters.pdb", ffOPTWR }
1106   };
1107 #define NFILE asize(fnm)
1108
1109   CopyRight(stderr,argv[0]);
1110   parse_common_args(&argc,argv,
1111                     PCA_CAN_VIEW | PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
1112                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,
1113                     &oenv);
1114
1115   /* parse options */
1116   bReadMat   = opt2bSet("-dm",NFILE,fnm);
1117   bReadTraj  = opt2bSet("-f",NFILE,fnm) || !bReadMat;
1118   if ( opt2parg_bSet("-av",asize(pa),pa) ||
1119        opt2parg_bSet("-wcl",asize(pa),pa) ||
1120        opt2parg_bSet("-nst",asize(pa),pa) ||
1121        opt2parg_bSet("-rmsmin",asize(pa),pa) ||
1122        opt2bSet("-cl",NFILE,fnm) )
1123     trx_out_fn = opt2fn("-cl",NFILE,fnm);
1124   else
1125     trx_out_fn = NULL;
1126   if (bReadMat && output_env_get_time_factor(oenv)!=1) {
1127     fprintf(stderr,
1128             "\nWarning: assuming the time unit in %s is %s\n",
1129             opt2fn("-dm",NFILE,fnm),output_env_get_time_unit(oenv));
1130   }
1131   if (trx_out_fn && !bReadTraj)
1132     fprintf(stderr,"\nWarning: "
1133             "cannot write cluster structures without reading trajectory\n"
1134             "         ignoring option -cl %s\n", trx_out_fn);
1135
1136   method=1;
1137   while ( method < m_nr && gmx_strcasecmp(methodname[0], methodname[method])!=0 )
1138     method++;
1139   if (method == m_nr)
1140     gmx_fatal(FARGS,"Invalid method");
1141
1142   bAnalyze = (method == m_linkage || method == m_jarvis_patrick ||
1143               method == m_gromos );
1144
1145   /* Open log file */
1146   log = ftp2FILE(efLOG,NFILE,fnm,"w");
1147
1148   fprintf(stderr,"Using %s method for clustering\n",methodname[0]);
1149   fprintf(log,"Using %s method for clustering\n",methodname[0]);
1150
1151   /* check input and write parameters to log file */
1152   bUseRmsdCut = FALSE;
1153   if (method == m_jarvis_patrick) {
1154     bJP_RMSD = (M == 0) || opt2parg_bSet("-cutoff",asize(pa),pa);
1155     if ((M<0) || (M == 1))
1156       gmx_fatal(FARGS,"M (%d) must be 0 or larger than 1",M);
1157     if (M < 2) {
1158       sprintf(buf1,"Will use P=%d and RMSD cutoff (%g)",P,rmsdcut);
1159       bUseRmsdCut = TRUE;
1160     } else {
1161       if (P >= M)
1162         gmx_fatal(FARGS,"Number of neighbors required (P) must be less than M");
1163       if (bJP_RMSD) {
1164         sprintf(buf1,"Will use P=%d, M=%d and RMSD cutoff (%g)",P,M,rmsdcut);
1165         bUseRmsdCut = TRUE;
1166       } else
1167         sprintf(buf1,"Will use P=%d, M=%d",P,M);
1168     }
1169     ffprintf1(stderr,log,buf,"%s for determining the neighbors\n\n",buf1);
1170   } else /* method != m_jarvis */
1171     bUseRmsdCut = ( bBinary || method == m_linkage || method == m_gromos );
1172   if (bUseRmsdCut && method != m_jarvis_patrick)
1173     fprintf(log,"Using RMSD cutoff %g nm\n",rmsdcut);
1174   if ( method==m_monte_carlo )
1175     fprintf(log,"Using %d iterations\n",niter);
1176
1177   if (skip < 1)
1178     gmx_fatal(FARGS,"skip (%d) should be >= 1",skip);
1179
1180   /* get input */
1181   if (bReadTraj) {
1182     /* don't read mass-database as masses (and top) are not used */
1183     read_tps_conf(ftp2fn(efTPS,NFILE,fnm),buf,&top,&ePBC,&xtps,NULL,box,
1184                   bAnalyze);
1185
1186     fprintf(stderr,"\nSelect group for least squares fit%s:\n",
1187             bReadMat?"":" and RMSD calculation");
1188     get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
1189               1,&ifsize,&fitidx,&grpname);
1190     if (trx_out_fn) {
1191       fprintf(stderr,"\nSelect group for output:\n");
1192       get_index(&(top.atoms),ftp2fn_null(efNDX,NFILE,fnm),
1193                 1,&iosize,&outidx,&grpname);
1194       /* merge and convert both index groups: */
1195       /* first copy outidx to index. let outidx refer to elements in index */
1196       snew(index,iosize);
1197       isize = iosize;
1198       for(i=0; i<iosize; i++) {
1199         index[i]=outidx[i];
1200         outidx[i]=i;
1201       }
1202       /* now lookup elements from fitidx in index, add them if necessary
1203          and also let fitidx refer to elements in index */
1204       for(i=0; i<ifsize; i++) {
1205         j=0;
1206         while (j<isize && index[j]!=fitidx[i])
1207           j++;
1208         if (j>=isize) {
1209           /* slow this way, but doesn't matter much */
1210           isize++;
1211           srenew(index,isize);
1212         }
1213         index[j]=fitidx[i];
1214         fitidx[i]=j;
1215       }
1216     } else { /* !trx_out_fn */
1217       isize = ifsize;
1218       snew(index, isize);
1219       for(i=0; i<ifsize; i++) {
1220         index[i]=fitidx[i];
1221         fitidx[i]=i;
1222       }
1223     }
1224   }
1225   /* Initiate arrays */
1226   snew(d1,isize);
1227   snew(d2,isize);
1228   for(i=0; (i<isize); i++) {
1229     snew(d1[i],isize);
1230     snew(d2[i],isize);
1231   }
1232
1233   if (bReadTraj) {
1234     /* Loop over first coordinate file */
1235     fn = opt2fn("-f",NFILE,fnm);
1236
1237     xx = read_whole_trj(fn,isize,index,skip,&nf,&time,oenv);
1238     output_env_conv_times(oenv, nf, time);
1239     if (!bRMSdist || bAnalyze) {
1240       /* Center all frames on zero */
1241       snew(mass,isize);
1242       for(i=0; i<ifsize; i++)
1243         mass[fitidx[i]] = top.atoms.atom[index[fitidx[i]]].m;
1244       if (bFit)
1245       for(i=0; i<nf; i++)
1246         reset_x(ifsize,fitidx,isize,NULL,xx[i],mass);
1247     }
1248   }
1249   if (bReadMat) {
1250     fprintf(stderr,"Reading rms distance matrix ");
1251     read_xpm_matrix(opt2fn("-dm",NFILE,fnm),&readmat);
1252     fprintf(stderr,"\n");
1253     if (readmat[0].nx != readmat[0].ny)
1254       gmx_fatal(FARGS,"Matrix (%dx%d) is not square",
1255                   readmat[0].nx,readmat[0].ny);
1256     if (bReadTraj && bAnalyze && (readmat[0].nx != nf))
1257       gmx_fatal(FARGS,"Matrix size (%dx%d) does not match the number of "
1258                   "frames (%d)",readmat[0].nx,readmat[0].ny,nf);
1259
1260     nf = readmat[0].nx;
1261     sfree(time);
1262     time = readmat[0].axis_x;
1263     time_invfac = output_env_get_time_invfactor(oenv);
1264     for(i=0; i<nf; i++)
1265       time[i] *= time_invfac;
1266
1267     rms = init_mat(readmat[0].nx,method == m_diagonalize);
1268     convert_mat(&(readmat[0]),rms);
1269
1270     nlevels = readmat[0].nmap;
1271   } else { /* !bReadMat */
1272     rms = init_mat(nf,method == m_diagonalize);
1273     nrms = (nf*(nf-1))/2;
1274     if (!bRMSdist) {
1275       fprintf(stderr,"Computing %dx%d RMS deviation matrix\n",nf,nf);
1276       snew(x1,isize);
1277       for(i1=0; (i1<nf); i1++) {
1278         for(i2=i1+1; (i2<nf); i2++) {
1279           for(i=0; i<isize; i++)
1280             copy_rvec(xx[i1][i],x1[i]);
1281           if (bFit)
1282             do_fit(isize,mass,xx[i2],x1);
1283           rmsd = rmsdev(isize,mass,xx[i2],x1);
1284           set_mat_entry(rms,i1,i2,rmsd);
1285         }
1286         nrms -= (nf-i1-1);
1287         fprintf(stderr,"\r# RMSD calculations left: %d   ",nrms);
1288       }
1289     } else { /* bRMSdist */
1290       fprintf(stderr,"Computing %dx%d RMS distance deviation matrix\n",nf,nf);
1291       for(i1=0; (i1<nf); i1++) {
1292         calc_dist(isize,xx[i1],d1);
1293         for(i2=i1+1; (i2<nf); i2++) {
1294           calc_dist(isize,xx[i2],d2);
1295           set_mat_entry(rms,i1,i2,rms_dist(isize,d1,d2));
1296       }
1297         nrms -= (nf-i1-1);
1298         fprintf(stderr,"\r# RMSD calculations left: %d   ",nrms);
1299       }
1300     }
1301     fprintf(stderr,"\n\n");
1302   }
1303   ffprintf2(stderr,log,buf,"The RMSD ranges from %g to %g nm\n",
1304             rms->minrms,rms->maxrms);
1305   ffprintf1(stderr,log,buf,"Average RMSD is %g\n",2*rms->sumrms/(nf*(nf-1)));
1306   ffprintf1(stderr,log,buf,"Number of structures for matrix %d\n",nf);
1307   ffprintf1(stderr,log,buf,"Energy of the matrix is %g nm\n",mat_energy(rms));
1308   if (bUseRmsdCut && (rmsdcut < rms->minrms || rmsdcut > rms->maxrms) )
1309     fprintf(stderr,"WARNING: rmsd cutoff %g is outside range of rmsd values "
1310             "%g to %g\n",rmsdcut,rms->minrms,rms->maxrms);
1311   if (bAnalyze && (rmsmin < rms->minrms) )
1312     fprintf(stderr,"WARNING: rmsd minimum %g is below lowest rmsd value %g\n",
1313             rmsmin,rms->minrms);
1314   if (bAnalyze && (rmsmin > rmsdcut) )
1315     fprintf(stderr,"WARNING: rmsd minimum %g is above rmsd cutoff %g\n",
1316             rmsmin,rmsdcut);
1317
1318   /* Plot the rmsd distribution */
1319   rmsd_distribution(opt2fn("-dist",NFILE,fnm),rms,oenv);
1320
1321   if (bBinary) {
1322     for(i1=0; (i1 < nf); i1++)
1323       for(i2=0; (i2 < nf); i2++)
1324         if (rms->mat[i1][i2] < rmsdcut)
1325           rms->mat[i1][i2] = 0;
1326         else
1327           rms->mat[i1][i2] = 1;
1328   }
1329
1330   snew(clust.cl,nf);
1331   switch (method) {
1332   case m_linkage:
1333     /* Now sort the matrix and write it out again */
1334     gather(rms,rmsdcut,&clust);
1335     break;
1336   case m_diagonalize:
1337     /* Do a diagonalization */
1338       snew(eigval,nf);
1339       snew(tmp,nf*nf);
1340       memcpy(tmp,rms->mat[0],nf*nf*sizeof(real));
1341       eigensolver(tmp,nf,0,nf,eigval,rms->mat[0]);
1342       sfree(tmp);
1343
1344       fp = xvgropen(opt2fn("-ev",NFILE,fnm),"RMSD matrix Eigenvalues",
1345                     "Eigenvector index","Eigenvalues (nm\\S2\\N)",oenv);
1346       for(i=0; (i<nf); i++)
1347           fprintf(fp,"%10d  %10g\n",i,eigval[i]);
1348           ffclose(fp);
1349       break;
1350   case m_monte_carlo:
1351     mc_optimize(log,rms,niter,&seed,kT);
1352     swap_mat(rms);
1353     reset_index(rms);
1354     break;
1355   case m_jarvis_patrick:
1356     jarvis_patrick(rms->nn,rms->mat,M,P,bJP_RMSD ? rmsdcut : -1,&clust);
1357     break;
1358   case m_gromos:
1359     gromos(rms->nn,rms->mat,rmsdcut,&clust);
1360     break;
1361   default:
1362     gmx_fatal(FARGS,"DEATH HORROR unknown method \"%s\"",methodname[0]);
1363   }
1364
1365   if (method == m_monte_carlo || method == m_diagonalize)
1366     fprintf(stderr,"Energy of the matrix after clustering is %g nm\n",
1367             mat_energy(rms));
1368
1369   if (bAnalyze) {
1370     if (minstruct > 1) {
1371       ncluster = plot_clusters(nf,rms->mat,&clust,nlevels,minstruct);
1372     } else {
1373       mark_clusters(nf,rms->mat,rms->maxrms,&clust);
1374     }
1375     init_t_atoms(&useatoms,isize,FALSE);
1376     snew(usextps, isize);
1377     useatoms.resinfo = top.atoms.resinfo;
1378     for(i=0; i<isize; i++) {
1379       useatoms.atomname[i]=top.atoms.atomname[index[i]];
1380       useatoms.atom[i].resind = top.atoms.atom[index[i]].resind;
1381       useatoms.nres = max(useatoms.nres,useatoms.atom[i].resind+1);
1382       copy_rvec(xtps[index[i]],usextps[i]);
1383     }
1384     useatoms.nr=isize;
1385     analyze_clusters(nf,&clust,rms->mat,isize,&useatoms,usextps,mass,xx,time,
1386                      ifsize,fitidx,iosize,outidx,
1387                      bReadTraj?trx_out_fn:NULL,
1388                      opt2fn_null("-sz",NFILE,fnm),
1389                      opt2fn_null("-tr",NFILE,fnm),
1390                      opt2fn_null("-ntr",NFILE,fnm),
1391                      opt2fn_null("-clid",NFILE,fnm),
1392                      bAverage, write_ncl, write_nst, rmsmin, bFit, log,
1393                      rlo_bot,rhi_bot,oenv);
1394   }
1395   ffclose(log);
1396
1397   if (bBinary && !bAnalyze)
1398     /* Make the clustering visible */
1399     for(i2=0; (i2 < nf); i2++)
1400        for(i1=i2+1; (i1 < nf); i1++)
1401          if (rms->mat[i1][i2])
1402            rms->mat[i1][i2] = rms->maxrms;
1403
1404   fp = opt2FILE("-o",NFILE,fnm,"w");
1405   fprintf(stderr,"Writing rms distance/clustering matrix ");
1406   if (bReadMat) {
1407     write_xpm(fp,0,readmat[0].title,readmat[0].legend,readmat[0].label_x,
1408               readmat[0].label_y,nf,nf,readmat[0].axis_x,readmat[0].axis_y,
1409               rms->mat,0.0,rms->maxrms,rlo_top,rhi_top,&nlevels);
1410   }
1411   else {
1412     sprintf(buf,"Time (%s)",output_env_get_time_unit(oenv));
1413     sprintf(title,"RMS%sDeviation / Cluster Index",
1414             bRMSdist ? " Distance " : " ");
1415     if (minstruct > 1) {
1416       write_xpm_split(fp,0,title,"RMSD (nm)",buf,buf,
1417                       nf,nf,time,time,rms->mat,0.0,rms->maxrms,&nlevels,
1418                       rlo_top,rhi_top,0.0,(real) ncluster,
1419                       &ncluster,TRUE,rlo_bot,rhi_bot);
1420     } else {
1421       write_xpm(fp,0,title,"RMSD (nm)",buf,buf,
1422                 nf,nf,time,time,rms->mat,0.0,rms->maxrms,
1423                 rlo_top,rhi_top,&nlevels);
1424     }
1425   }
1426   fprintf(stderr,"\n");
1427   ffclose(fp);
1428
1429   /* now show what we've done */
1430   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
1431   do_view(oenv,opt2fn_null("-sz",NFILE,fnm),"-nxy");
1432   if (method == m_diagonalize)
1433     do_view(oenv,opt2fn_null("-ev",NFILE,fnm),"-nxy");
1434   do_view(oenv,opt2fn("-dist",NFILE,fnm),"-nxy");
1435   if (bAnalyze) {
1436     do_view(oenv,opt2fn_null("-tr",NFILE,fnm),"-nxy");
1437     do_view(oenv,opt2fn_null("-ntr",NFILE,fnm),"-nxy");
1438     do_view(oenv,opt2fn_null("-clid",NFILE,fnm),"-nxy");
1439   }
1440
1441   /* Thank the user for her patience */
1442   thanx(stderr);
1443
1444   return 0;
1445 }