Project

General

Profile

0001-Add-new-tool-g_correlation.-This-tool-computes-Pears.patch

patch againts git master - Alexey Shvetsov, 06/15/2010 01:31 AM

View differences:

src/tools/CMakeLists.txt
24 24
            gmx_sgangle.c   gmx_sorient.c   gmx_spol.c      gmx_tcaf.c      
25 25
            gmx_traj.c      gmx_velacc.c    gmx_helixorient.c 
26 26
            gmx_clustsize.c gmx_mdmat.c     gmx_wham.c      
27
            correl.c        gmx_sham.c      gmx_nmtraj.c    
27
            correl.c        gmx_sham.c      gmx_nmtraj.c    gmx_correlation.c
28 28
            gmx_trjconv.c   gmx_trjcat.c    gmx_trjorder.c  gmx_xpm2ps.c    
29 29
            gmx_editconf.c  gmx_genbox.c    gmx_genion.c    gmx_genconf.c   
30 30
            gmx_genpr.c     gmx_eneconv.c   gmx_vanhove.c   gmx_wheel.c     
......
43 43
    g_angle g_bond g_bundle g_chi g_cluster g_confrms g_covar
44 44
    g_current g_density g_densmap g_dih g_dielectric
45 45
    g_helixorient g_principal g_dipoles g_disre g_dist
46
    g_dyndom g_enemat g_energy g_lie g_filter g_gyrate
46
    g_dyndom g_enemat g_energy g_lie g_filter g_gyrate g_correlation
47 47
    g_h2order g_hbond g_helix g_mindist g_msd g_morph g_nmeig
48 48
    g_nmens g_order g_polystat g_potential g_rama g_rdf g_rms
49 49
    g_rmsf g_rotacf g_saltbr g_sas g_select g_sgangle g_sham g_sorient
src/tools/Makefile.am
44 44
	gmx_editconf.c	gmx_genbox.c	gmx_genion.c	gmx_genconf.c	\
45 45
	gmx_genpr.c	gmx_eneconv.c	gmx_vanhove.c	gmx_wheel.c	\
46 46
	addconf.c 	addconf.h	gmx_tune_pme.c    \
47
	calcpot.c 	calcpot.h 	edittop.c
47
	calcpot.c 	calcpot.h 	edittop.c	gmx_correlation.c
48 48

  
49 49
bin_PROGRAMS = \
50 50
	do_dssp		editconf	eneconv		\
src/tools/g_correlation.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 <gmx_ana.h>
40

  
41

  
42

  
43
/* This is just a wrapper binary.
44
* The code that used to be in g_covar.c is now in gmx_covar.c,
45
* where the old main function is called gmx_covar().
46
*/
47
int
48
main(int argc, char *argv[])
49
{
50
  gmx_correlation(argc,argv);
51
  return 0;
52
}
53

  
54

  
55
  
src/tools/gmx_correlation.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
#include <math.h>
39
#include <string.h>
40
#include <time.h>
41

  
42
#if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
43
#include <direct.h>
44
#include <io.h>
45
#endif
46

  
47
#include "statutil.h"
48
#include "sysstuff.h"
49
#include "typedefs.h"
50
#include "smalloc.h"
51
#include "macros.h"
52
#include "vec.h"
53
#include "pbc.h"
54
#include "copyrite.h"
55
#include "futil.h"
56
#include "statutil.h"
57
#include "index.h"
58
#include "confio.h"
59
#include "trnio.h"
60
#include "mshift.h"
61
#include "xvgr.h"
62
#include "do_fit.h"
63
#include "rmpbc.h"
64
#include "txtdump.h"
65
#include "matio.h"
66
#include "physics.h"
67
#include "gmx_ana.h"
68

  
69
int gmx_correlation(int argc,char *argv[])
70
{
71
  const char *desc[] = {
72
    "[TT]g_correlation[tt] calculates and the correlation matrix.",
73
    "All structures are fitted to the structure in the structure file.",
74
    "When this is not a run input file periodicity will not be taken into",
75
    "account.",
76
    "[PAR]",
77
    "When the same atoms are used for the fit and the correlation analysis,",
78
    "the reference structure for the fit is written first with t=-1.",
79
    "[PAR]",
80
    "Option [TT]-ascii[tt] writes the whole correlation matrix to",
81
    "an ASCII file.", 
82
    "[PAR]",
83
    "Option [TT]-xpm[tt] writes the whole covariance matrix to an xpm file.",
84
    "[PAR]"
85
  };
86
  static bool bFit=TRUE,bRef=FALSE,bPBC=TRUE;
87
  t_pargs pa[] = {
88
    { "-fit",  FALSE, etBOOL, {&bFit},
89
      "Fit to a reference structure"},
90
    { "-ref",  FALSE, etBOOL, {&bRef},
91
      "Use the deviation from the conformation in the structure file instead of from the average" },
92
    { "-pbc",  FALSE,  etBOOL, {&bPBC},
93
      "Apply corrections for periodic boundary conditions" }
94
  };
95
  FILE       *out;
96
  int        status;
97
  t_topology top;
98
  int        ePBC;
99
  t_atoms    *atoms;  
100
  rvec       *xread,*xref;
101
  rvec       *xa,*xava;
102
  rvec       *xb,*xavb;
103
  matrix     box,zerobox;
104
  real       inv_nframes;
105
  real       t,tstart,tend;
106
  real       *w_rls=NULL;
107
  real       min,max,*axisa,*axisb;
108
  real       **mat;
109
  real       *r2a,*r2b;
110
  int        natoms,nat,nframes0,nframes,nlevels;
111
  int        natomsa,natomsb;
112
  gmx_large_int_t i,j;
113
  const char *fitfile,*trxfile,*ndxfile;
114
  const char *asciifile,*xpmfile,*averfilea,*averfileb;
115
  char       str[STRLEN],*fitname;
116
  char       *ananamea,*ananameb;
117
  int        d,nfit;
118
  atom_id    *ifit;
119
  atom_id    *indexa,*indexb;
120
  t_rgb      rlo,rmi,rhi;
121
  output_env_t oenv;
122

  
123
  t_filenm fnm[] = { 
124
    { efTRX, "-f",  NULL, ffREAD }, 
125
    { efTPS, NULL,  NULL, ffREAD },
126
    { efNDX, NULL,  NULL, ffOPTRD },
127
    { efSTO, "-ava", "averagea.pdb", ffWRITE },
128
    { efSTO, "-avb", "averageb.pdb", ffWRITE },
129
    { efDAT, "-ascii","correlation", ffOPTWR },
130
    { efXPM, "-xpm","correlation", ffOPTWR }
131
  }; 
132
#define NFILE asize(fnm) 
133

  
134
  CopyRight(stderr,argv[0]); 
135
  parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
136
		    NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv); 
137

  
138
  clear_mat(zerobox);
139

  
140
  fitfile    = ftp2fn(efTPS,NFILE,fnm);
141
  trxfile    = ftp2fn(efTRX,NFILE,fnm);
142
  ndxfile    = ftp2fn_null(efNDX,NFILE,fnm);
143
  averfilea  = ftp2fn(efSTO,NFILE,fnm);
144
  averfileb  = ftp2fn(efSTO,NFILE,fnm);
145
  asciifile  = opt2fn_null("-ascii",NFILE,fnm);
146
  xpmfile    = opt2fn_null("-xpm",NFILE,fnm);
147

  
148
  read_tps_conf(fitfile,str,&top,&ePBC,&xref,NULL,box,TRUE);
149
  atoms=&top.atoms;
150

  
151
  if (bFit) {
152
    printf("\nChoose a group for the least squares fit\n"); 
153
    get_index(atoms,ndxfile,1,&nfit,&ifit,&fitname);
154
    if (nfit < 3) 
155
      gmx_fatal(FARGS,"Need >= 3 points to fit!\n");
156
  } else
157
    nfit=0;
158
  printf("\nChoose two groups for the correlation analysis\n"); 
159
  get_index(atoms,ndxfile,1,&natomsa,&indexa,&ananamea);
160
  get_index(atoms,ndxfile,1,&natomsb,&indexb,&ananameb);
161

  
162
  if (bFit) {
163
    snew(w_rls,atoms->nr);
164
    for(i=0; (i<nfit); i++)
165
	w_rls[ifit[i]]=1.0;
166
    }
167

  
168
  /* Prepare reference frame */
169
  if (bPBC)
170
   rm_pbc(&(top.idef),ePBC,atoms->nr,box,xref,xref);
171
  if (bFit)
172
   reset_x(nfit,ifit,atoms->nr,NULL,xref,w_rls);
173

  
174
  snew(xa,natomsa);
175
  snew(xb,natomsb);
176
  snew(xava,natomsa);
177
  snew(xavb,natomsb);
178
/*  if (sqrt(LARGE_INT_MAX)<natoms) {
179
    gmx_fatal(FARGS,"Number of degrees of freedoms to large for matrix.\n");
180
  } */
181

  
182
  fprintf(stderr,"Calculating the average structure ...\n");
183
  nframes0 = 0;
184
  nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
185
  if (nat != atoms->nr)
186
    fprintf(stderr,"\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n",natoms,nat);
187
  do {
188
    nframes0++;
189
    /* calculate x: a fitted struture of the selected atoms */
190
    if (bPBC)
191
      rm_pbc(&(top.idef),ePBC,nat,box,xread,xread);
192
    if (bFit) {
193
      reset_x(nfit,ifit,nat,NULL,xread,w_rls);
194
      do_fit(nat,w_rls,xref,xread);
195
    }
196
    for (i=0; i<natomsa; i++)
197
      rvec_inc(xava[i],xread[indexa[i]]);
198
    for (i=0; i<natomsb; i++)
199
      rvec_inc(xavb[i],xread[indexb[i]]);
200
  } while (read_next_x(oenv,status,&t,nat,xread,box));
201
  close_trj(status);
202
  
203
  inv_nframes = 1.0/nframes0;
204
  for(i=0; i<natomsa; i++)
205
    for(d=0; d<DIM; d++) {
206
      xava[i][d] *= inv_nframes;
207
      xread[indexa[i]][d] = xava[i][d];
208
    }
209
  write_sto_conf_indexed(opt2fn("-ava",NFILE,fnm),"Average structure",
210
  		         atoms,xread,NULL,epbcNONE,zerobox,natomsa,indexa);
211
  for(i=0; i<natomsb; i++)
212
   for(d=0; d<DIM; d++) {
213
      xavb[i][d] *= inv_nframes;
214
      xread[indexb[i]][d] = xavb[i][d];
215
   }
216
  write_sto_conf_indexed(opt2fn("-avb",NFILE,fnm),"Average structure",
217
			 atoms,xread,NULL,epbcNONE,zerobox,natomsb,indexb);
218
  sfree(xread);
219

  
220
  fprintf(stderr,"Constructing correlation matrix (%dx%d) ...\n",(int)natomsa,(int)natomsb);
221
  nframes=0;
222
  nat=read_first_x(oenv,&status,trxfile,&t,&xread,box);
223
  tstart = t;
224
  snew(r2a,natomsa);
225
  snew(r2b,natomsb);
226
  /* Allocate 2d matrix */
227
  snew(mat,natomsa);
228
  for (i=0;i<natomsb;i++)
229
  	snew(mat[i],natomsb);
230
  do {
231
    nframes++;
232
    tend = t;
233
    /* calculate x: a (fitted) structure of the selected atoms */
234
    if (bPBC)
235
      rm_pbc(&(top.idef),ePBC,nat,box,xread,xread);
236
    if (bFit) {
237
      reset_x(nfit,ifit,nat,NULL,xread,w_rls);
238
      do_fit(nat,w_rls,xref,xread);
239
    }
240
    if (bRef) {
241
      for (i=0; i<natomsa; i++)
242
	rvec_sub(xread[indexa[i]],xref[indexa[i]],xa[i]);
243
      for (i=0; i<natomsb;i++)
244
        rvec_sub(xread[indexb[i]],xref[indexb[i]],xb[i]);
245
      }
246
    else {
247
      for (i=0; i<natomsa; i++)
248
	rvec_sub(xread[indexa[i]],xava[i],xa[i]);
249
      for (i=0; i<natomsb; i++)
250
        rvec_sub(xread[indexb[i]],xavb[i],xb[i]);
251
      }
252
    /* find actual rmsf */
253
    for (i=0; i<natomsa; i++)
254
        r2a[i]+=norm2(xa[i]);
255
    for (i=0; i<natomsb; i++)
256
        r2b[i]+=norm2(xb[i]);
257
    /* TODO: Optimize calculation of correlation matrix */
258
    for (i=0; i<natomsa; i++) {
259
    	for (j=0; j<natomsb; j++) {
260
		mat[i][j] += iprod(xa[i],xb[j]);
261
	}
262
    }
263
  } while (read_next_x(oenv,status,&t,nat,xread,box) && 
264
	   (bRef || nframes < nframes0));
265
  close_trj(status);
266
  fprintf(stderr,"Read %d frames\n",nframes);
267
  /* Correct mat so it will be real correlation one */
268
  for (i=0;i<natomsa;i++) {
269
  	for (j=0;j<natomsb;j++) {
270
		mat[i][j] = mat[i][j]/sqrt(r2a[i]*r2b[j]);
271
	}
272
  }
273
  if (asciifile) {
274
    out = ffopen(asciifile,"w");
275
    for (i=0; i<natomsa; i++) {
276
      for (j=0; j<natomsb; j++) {
277
	fprintf(out,"%g ",mat[i][j]);
278
	}
279
      fprintf(out,"\n");
280
    }
281
    ffclose(out);
282
  }
283
  
284
  if (xpmfile) {
285
    min = 0;
286
    max = 0;
287
    for (i=0; i<natomsa; i++) {
288
      for (j=0; j<natomsb; j++) {
289
	if (mat[i][j] < min)
290
	  min = mat[i][j];
291
	if (mat[i][j] > max)
292
	  max = mat[i][j];
293
      }
294
    }
295
    snew(axisa,natomsa);
296
    snew(axisb,natomsb);
297
    for(i=0;i<natomsa;i++)
298
      axisa[i] = i+1;
299
    for(i=0;i<natomsb;i++)
300
      axisb[i] = i+1;
301
    rlo.r = 0; rlo.g = 0; rlo.b = 1;
302
    rmi.r = 1; rmi.g = 1; rmi.b = 1;
303
    rhi.r = 1; rhi.g = 0; rhi.b = 0;
304
    out = ffopen(xpmfile,"w");
305
    nlevels = 80;
306
    write_xpm3(out,0,"Correlation", "",
307
	       "Residue Number","Residue Number",natomsa,natomsb,axisa,axisb,
308
	       mat,min,0.0,max,rlo,rmi,rhi,&nlevels);
309
    ffclose(out);
310
    sfree(axisa);
311
    sfree(axisb);
312
  }
313
  for (i=0;i<natomsb;i++)
314
  	sfree(mat[i]);
315
  sfree(mat);
316
  sfree(r2a);
317
  sfree(r2b);
318

  
319
  thanx(stderr);
320

  
321
  return 0;
322
}
0
-