src/gromacs/gmxana/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, gmx_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     {
  79         srenew(m->mat[i], m->nn+deltan);
  80         for (j = m->nn; (j < m->nn+deltan); j++)
  81         {
  82             m->mat[i][j] = 0;
  83         }
  84     }
  85     /* Allocate new rows of the matrix, set energies to zero */
  86     for (i = m->nn; (i < m->nn+deltan); i++)
  87     {
  88         m->erow[i]  = 0;
  89         snew(m->mat[i], m->nn+deltan);
  90     }
  91     m->nn += deltan;
  92 }
  93
  94 void reset_index(t_mat *m)
  95 {
  96     int i;
  97
  98     for (i = 0; (i < m->n1); i++)
  99     {
 100         m->m_ind[i] = i;
 101     }
 102 }
 103
 104 void set_mat_entry(t_mat *m, int i, int j, real val)
 105 {
 106     m->mat[i][j] = m->mat[j][i] = val;
 107     m->maxrms    = max(m->maxrms, val);
 108     if (j != i)
 109     {
 110         m->minrms  = min(m->minrms, val);
 111     }
 112     m->sumrms   += val;
 113     m->nn        = max(m->nn, max(j+1, i+1));
 114 }
 115
 116 void done_mat(t_mat **m)
 117 {
 118     done_matrix((*m)->n1, &((*m)->mat));
 119     sfree((*m)->m_ind);
 120     sfree((*m)->erow);
 121     sfree(*m);
 122     *m = NULL;
 123 }
 124
 125 real row_energy(int nn, int row, real *mat)
 126 {
 127     real re = 0;
 128     int  i;
 129
 130     for (i = 0; (i < nn); i++)
 131     {
 132         re += abs(i-row)*mat[i];
 133     }
 134     return re/nn;
 135 }
 136
 137 real mat_energy(t_mat *m)
 138 {
 139     real re, retot;
 140     int  j, jj;
 141
 142     retot = 0;
 143     for (j = 0; (j < m->nn); j++)
 144     {
 145         jj         = m->m_ind[j];
 146         re         = row_energy(m->nn, jj, m->mat[j]);
 147         m->erow[j] = re;
 148         retot     += re;
 149     }
 150     m->emat = retot/m->nn;
 151     return m->emat;
 152 }
 153
 154 void swap_rows(t_mat *m, int isw, int jsw)
 155 {
 156     real *tmp, ttt;
 157     int   i;
 158
 159     /* Swap rows */
 160     tmp         = m->mat[isw];
 161     m->mat[isw] = m->mat[jsw];
 162     m->mat[jsw] = tmp;
 163     /* Swap columns */
 164     for (i = 0; (i < m->nn); i++)
 165     {
 166         ttt            = m->mat[isw][i];
 167         m->mat[isw][i] = m->mat[jsw][i];
 168         m->mat[jsw][i] = ttt;
 169     }
 170 }
 171
 172 void swap_mat(t_mat *m)
 173 {
 174     t_mat *tmp;
 175     int    i, j;
 176
 177     tmp = init_mat(m->nn, FALSE);
 178     for (i = 0; (i < m->nn); i++)
 179     {
 180         for (j = 0; (j < m->nn); j++)
 181         {
 182             tmp->mat[m->m_ind[i]][m->m_ind[j]] = m->mat[i][j];
 183         }
 184     }
 185     /*tmp->mat[i][j] =  m->mat[m->m_ind[i]][m->m_ind[j]]; */
 186     for (i = 0; (i < m->nn); i++)
 187     {
 188         for (j = 0; (j < m->nn); j++)
 189         {
 190             m->mat[i][j] = tmp->mat[i][j];
 191         }
 192     }
 193     done_mat(&tmp);
 194 }
 195
 196 void low_rmsd_dist(const char *fn, real maxrms, int nn, real **mat,
 197                    const output_env_t oenv)
 198 {
 199     FILE   *fp;
 200     int     i, j, *histo, x;
 201     real    fac;
 202
 203     fac = 100/maxrms;
 204     snew(histo, 101);
 205     for (i = 0; i < nn; i++)
 206     {
 207         for (j = i+1; j < nn; j++)
 208         {
 209             x = (int)(fac*mat[i][j]+0.5);
 210             if (x <= 100)
 211             {
 212                 histo[x]++;
 213             }
 214         }
 215     }
 216
 217     fp = xvgropen(fn, "RMS Distribution", "RMS (nm)", "a.u.", oenv);
 218     for (i = 0; (i < 101); i++)
 219     {
 220         fprintf(fp, "%10g  %10d\n", i/fac, histo[i]);
 221     }
 222     ffclose(fp);
 223     sfree(histo);
 224 }
 225
 226 void rmsd_distribution(const char *fn, t_mat *rms, const output_env_t oenv)
 227 {
 228     low_rmsd_dist(fn, rms->maxrms, rms->nn, rms->mat, oenv);
 229 }
 230
 231 t_clustid *new_clustid(int n1)
 232 {
 233     t_clustid *c;
 234     int        i;
 235
 236     snew(c, n1);
 237     for (i = 0; (i < n1); i++)
 238     {
 239         c[i].conf  = i;
 240         c[i].clust = i;
 241     }
 242     return c;
 243 }