src/tools/cmat.c

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