Fixed make_edi.c
[gromacs/rigid-bodies.git] / src / tools / cmat.c
blob4685c7c9990585db8d8d9a1421ddee91c72f4c62
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
39 #include "cmat.h"
40 #include "smalloc.h"
41 #include "macros.h"
42 #include "xvgr.h"
43 #include "matio.h"
44 #include "futil.h"
46 t_mat *init_mat(int n1,bool b1D)
48 t_mat *m;
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);
61 snew(m->erow,n1);
62 snew(m->m_ind,n1);
63 reset_index(m);
65 return m;
68 void enlarge_mat(t_mat *m,int deltan)
70 int i,j;
72 srenew(m->erow,m->nn+deltan);
73 srenew(m->m_ind,m->nn+deltan);
74 srenew(m->mat,m->nn+deltan);
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;
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);
87 m->nn += deltan;
90 void reset_index(t_mat *m)
92 int i;
94 for(i=0; (i<m->n1); i++)
95 m->m_ind[i] = i;
98 void set_mat_entry(t_mat *m,int i,int j,real val)
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));
108 void done_mat(t_mat **m)
110 done_matrix((*m)->n1,&((*m)->mat));
111 sfree((*m)->m_ind);
112 sfree((*m)->erow);
113 sfree(*m);
114 *m = NULL;
117 real row_energy(int nn,int row,real *mat)
119 real re = 0;
120 int i;
122 for(i=0; (i<nn); i++) {
123 re += abs(i-row)*mat[i];
125 return re/nn;
128 real mat_energy(t_mat *m)
130 real re,retot;
131 int j,jj;
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;
140 m->emat = retot/m->nn;
141 return m->emat;
144 void swap_rows(t_mat *m,int isw,int jsw)
146 real *tmp,ttt;
147 int i;
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;
161 void swap_mat(t_mat *m)
163 t_mat *tmp;
164 int i,j;
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);
177 void low_rmsd_dist(const char *fn,real maxrms,int nn,real **mat,
178 const output_env_t oenv)
180 FILE *fp;
181 int i,j,*histo,x;
182 real fac;
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]++;
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);
200 void rmsd_distribution(const char *fn,t_mat *rms,const output_env_t oenv)
202 low_rmsd_dist(fn,rms->maxrms,rms->nn,rms->mat,oenv);
205 t_clustid *new_clustid(int n1)
207 t_clustid *c;
208 int i;
210 snew(c,n1);
211 for(i=0; (i<n1); i++) {
212 c[i].conf = i;
213 c[i].clust = i;
215 return c;