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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gromacs Runs One Microsecond At Cannonball Speeds
  31  */
  32
  33 #include "cmat.h"
  34 #include "smalloc.h"
  35 #include "macros.h"
  36 #include "xvgr.h"
  37 #include "matio.h"
  38
  39 t_mat *init_mat(int n1,bool b1D)
  40 {
  41   t_mat *m;
  42
  43   snew(m,1);
  44   m->n1   = n1;
  45   m->nn   = 0;
  46   m->b1D  = b1D;
  47   m->emat = 0;
  48   m->maxrms = 0;
  49   m->minrms = 1e20;
  50   m->sumrms = 0;
  51   m->nn   = 0;
  52   m->mat  = mk_matrix(n1,n1,b1D);
  53
  54   snew(m->erow,n1);
  55   snew(m->m_ind,n1);
  56   reset_index(m);
  57
  58   return m;
  59 }
  60
  61 void enlarge_mat(t_mat *m,int deltan)
  62 {
  63   int i,j;
  64
  65   srenew(m->erow,m->nn+deltan);
  66   srenew(m->m_ind,m->nn+deltan);
  67   srenew(m->mat,m->nn+deltan);
  68
  69   /* Reallocate existing rows in the matrix, and set them to zero */
  70   for(i=0; (i<m->nn); i++) {
  71     srenew(m->mat[i],m->nn+deltan);
  72     for(j=m->nn; (j<m->nn+deltan); j++)
  73       m->mat[i][j] = 0;
  74   }
  75   /* Allocate new rows of the matrix, set energies to zero */
  76   for(i=m->nn; (i<m->nn+deltan); i++) {
  77     m->erow[i]  = 0;
  78     snew(m->mat[i],m->nn+deltan);
  79   }
  80   m->nn += deltan;
  81 }
  82
  83 void reset_index(t_mat *m)
  84 {
  85   int i;
  86
  87   for(i=0; (i<m->n1); i++)
  88     m->m_ind[i] = i;
  89 }
  90
  91 void set_mat_entry(t_mat *m,int i,int j,real val)
  92 {
  93   m->mat[i][j] = m->mat[j][i] = val;
  94   m->maxrms    = max(m->maxrms,val);
  95   if (j!=i)
  96     m->minrms  = min(m->minrms,val);
  97   m->sumrms   += val;
  98   m->nn        = max(m->nn,max(j+1,i+1));
  99 }
 100
 101 void done_mat(t_mat **m)
 102 {
 103   done_matrix((*m)->n1,&((*m)->mat));
 104   sfree((*m)->m_ind);
 105   sfree((*m)->erow);
 106   sfree(*m);
 107   *m = NULL;
 108 }
 109
 110 real row_energy(int nn,int row,real *mat)
 111 {
 112   real re = 0;
 113   int  i;
 114
 115   for(i=0; (i<nn); i++) {
 116     re += abs(i-row)*mat[i];
 117   }
 118   return re/nn;
 119 }
 120
 121 real mat_energy(t_mat *m)
 122 {
 123   real re,retot;
 124   int  j,jj;
 125
 126   retot = 0;
 127   for(j=0; (j<m->nn); j++) {
 128     jj = m->m_ind[j];
 129     re = row_energy(m->nn,jj,m->mat[j]);
 130     m->erow[j] = re;
 131     retot += re;
 132   }
 133   m->emat = retot/m->nn;
 134   return m->emat;
 135 }
 136
 137 void swap_rows(t_mat *m,int isw,int jsw)
 138 {
 139   real *tmp,ttt;
 140   int  i;
 141
 142   /* Swap rows */
 143   tmp         = m->mat[isw];
 144   m->mat[isw] = m->mat[jsw];
 145   m->mat[jsw] = tmp;
 146   /* Swap columns */
 147   for(i=0; (i<m->nn); i++) {
 148     ttt            = m->mat[isw][i];
 149     m->mat[isw][i] = m->mat[jsw][i];
 150     m->mat[jsw][i] = ttt;
 151   }
 152 }
 153
 154 void swap_mat(t_mat *m)
 155 {
 156   t_mat *tmp;
 157   int   i,j;
 158
 159   tmp = init_mat(m->nn,FALSE);
 160   for(i=0; (i<m->nn); i++)
 161     for(j=0; (j<m->nn); j++)
 162       tmp->mat[m->m_ind[i]][m->m_ind[j]] = m->mat[i][j];
 163   /*tmp->mat[i][j] =  m->mat[m->m_ind[i]][m->m_ind[j]]; */
 164   for(i=0; (i<m->nn); i++)
 165     for(j=0; (j<m->nn); j++)
 166       m->mat[i][j] = tmp->mat[i][j];
 167   done_mat(&tmp);
 168 }
 169
 170 void low_rmsd_dist(char *fn,real maxrms,int nn,real **mat)
 171 {
 172   FILE   *fp;
 173   int    i,j,*histo,x;
 174   real   fac;
 175
 176   fac = 100/maxrms;
 177   snew(histo,101);
 178   for(i=0; i<nn; i++)
 179     for(j=i+1; j<nn; j++) {
 180       x = (int)(fac*mat[i][j]+0.5);
 181       if (x <= 100)
 182         histo[x]++;
 183     }
 184
 185   fp = xvgropen(fn,"RMS Distribution","RMS (nm)","a.u.");
 186   for(i=0; (i<101); i++)
 187     fprintf(fp,"%10g  %10d\n",i/fac,histo[i]);
 188   fclose(fp);
 189   sfree(histo);
 190 }
 191
 192 void rmsd_distribution(char *fn,t_mat *rms)
 193 {
 194   low_rmsd_dist(fn,rms->maxrms,rms->nn,rms->mat);
 195 }
 196
 197 t_clustid *new_clustid(int n1)
 198 {
 199   t_clustid *c;
 200   int       i;
 201
 202   snew(c,n1);
 203   for(i=0; (i<n1); i++) {
 204     c[i].conf = i;
 205     c[i].clust = i;
 206   }
 207   return c;
 208 }
 209