src/tools/cmat.c

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