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