src/tools/cmat.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_cmat_c = "$Id$";
  30
  31 #include "cmat.h"
  32 #include "smalloc.h"
  33 #include "macros.h"
  34 #include "xvgr.h"
  35
  36 real **mk_matrix(int n1,bool b1D)
  37 {
  38   real **m;
  39   int  i;
  40
  41   snew(m,n1);
  42   if (b1D)
  43     snew(m[0],n1*n1);
  44
  45   for(i=0; (i<n1); i++) {
  46     if (b1D)
  47       m[i] = &(m[0][i*n1]);
  48     else
  49       snew(m[i],n1);
  50   }
  51   return m;
  52 }
  53
  54 void done_matrix(int n1,real ***m)
  55 {
  56   int i;
  57
  58   for(i=0; (i<n1); i++)
  59     sfree((*m)[i]);
  60   sfree(*m);
  61   *m = NULL;
  62 }
  63
  64 t_mat *init_mat(int n1,bool b1D)
  65 {
  66   t_mat *m;
  67   int   i;
  68
  69   snew(m,1);
  70   m->n1   = n1;
  71   m->nn   = 0;
  72   m->b1D  = b1D;
  73   m->emat = 0;
  74   m->maxrms = 0;
  75   m->sumrms = 0;
  76   m->nn   = 0;
  77   m->mat  = mk_matrix(n1,b1D);
  78
  79   snew(m->erow,n1);
  80   snew(m->m_ind,n1);
  81   reset_index(m);
  82
  83   return m;
  84 }
  85
  86 void enlarge_mat(t_mat *m,int deltan)
  87 {
  88   int i,j;
  89
  90   srenew(m->erow,m->nn+deltan);
  91   srenew(m->m_ind,m->nn+deltan);
  92   srenew(m->mat,m->nn+deltan);
  93
  94   /* Reallocate existing rows in the matrix, and set them to zero */
  95   for(i=0; (i<m->nn); i++) {
  96     srenew(m->mat[i],m->nn+deltan);
  97     for(j=m->nn; (j<m->nn+deltan); j++)
  98       m->mat[i][j] = 0;
  99   }
 100   /* Allocate new rows of the matrix, set energies to zero */
 101   for(i=m->nn; (i<m->nn+deltan); i++) {
 102     m->erow[i]  = 0;
 103     snew(m->mat[i],m->nn+deltan);
 104   }
 105   m->nn += deltan;
 106 }
 107
 108 void reset_index(t_mat *m)
 109 {
 110   int i;
 111
 112   for(i=0; (i<m->n1); i++)
 113     m->m_ind[i] = i;
 114 }
 115
 116 void set_mat_entry(t_mat *m,int i,int j,real val)
 117 {
 118   m->mat[i][j] = m->mat[j][i] = val;
 119   m->maxrms    = max(m->maxrms,val);
 120   m->sumrms   += val;
 121   m->nn        = max(m->nn,max(j+1,i+1));
 122 }
 123
 124 void done_mat(t_mat **m)
 125 {
 126   int i;
 127
 128   done_matrix((*m)->n1,&((*m)->mat));
 129   sfree((*m)->m_ind);
 130   sfree((*m)->erow);
 131   sfree(*m);
 132   *m = NULL;
 133 }
 134
 135 real row_energy(int nn,int row,real *mat,int m_ind[])
 136 {
 137   real re = 0;
 138   int  i;
 139
 140   for(i=0; (i<nn); i++) {
 141     re += abs(i-row)*mat[i];
 142   }
 143   return re/nn;
 144 }
 145
 146 real mat_energy(t_mat *m)
 147 {
 148   real re,retot;
 149   int  j,jj;
 150
 151   retot = 0;
 152   for(j=0; (j<m->nn); j++) {
 153     jj = m->m_ind[j];
 154     re = row_energy(m->nn,jj,m->mat[j],m->m_ind);
 155     m->erow[j] = re;
 156     retot += re;
 157   }
 158   m->emat = retot/m->nn;
 159   return m->emat;
 160 }
 161
 162 void swap_rows(t_mat *m,int isw,int jsw)
 163 {
 164   real *tmp,ttt;
 165   int  i;
 166
 167   /* Swap rows */
 168   tmp         = m->mat[isw];
 169   m->mat[isw] = m->mat[jsw];
 170   m->mat[jsw] = tmp;
 171   /* Swap columns */
 172   for(i=0; (i<m->nn); i++) {
 173     ttt            = m->mat[isw][i];
 174     m->mat[isw][i] = m->mat[jsw][i];
 175     m->mat[jsw][i] = ttt;
 176   }
 177 }
 178
 179 void swap_mat(t_mat *m)
 180 {
 181   t_mat *tmp;
 182   int   i,j;
 183
 184   tmp = init_mat(m->nn,FALSE);
 185   for(i=0; (i<m->nn); i++)
 186     for(j=0; (j<m->nn); j++)
 187       tmp->mat[m->m_ind[i]][m->m_ind[j]] = m->mat[i][j];
 188   /*tmp->mat[i][j] =  m->mat[m->m_ind[i]][m->m_ind[j]]; */
 189   for(i=0; (i<m->nn); i++)
 190     for(j=0; (j<m->nn); j++)
 191       m->mat[i][j] = tmp->mat[i][j];
 192   done_mat(&tmp);
 193 }
 194
 195 void low_rmsd_dist(char *fn,real maxrms,int nn,real **mat)
 196 {
 197   FILE   *fp;
 198   int    i,j,*histo;
 199   real   fac;
 200
 201   fac = 100/maxrms;
 202   snew(histo,101);
 203   for(i=0; (i<nn); i++)
 204     for(j=i+1; (j<nn); j++)
 205       histo[(int)(fac*mat[i][j])]++;
 206
 207   fp = xvgropen(fn,"RMS Distribution","RMS (nm)","a.u.");
 208   for(i=0; (i<101); i++)
 209     fprintf(fp,"%10g  %10d\n",i/fac,histo[i]);
 210   fclose(fp);
 211   sfree(histo);
 212   xvgr_file(fn,NULL);
 213 }
 214
 215 void rmsd_distribution(char *fn,t_mat *rms)
 216 {
 217   low_rmsd_dist(fn,rms->maxrms,rms->nn,rms->mat);
 218 }
 219
 220 t_clustid *new_clustid(int n1)
 221 {
 222   t_clustid *c;
 223   int       i;
 224
 225   snew(c,n1);
 226   for(i=0; (i<n1); i++) {
 227     c[i].conf = i;
 228     c[i].clust = i;
 229   }
 230   return c;
 231 }
 232