Project

General

Profile

0001-Add-new-tool-g_correl.patch

new patch againts git master - Alexey Shvetsov, 06/15/2010 04:49 PM

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_correl.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_correl
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_correl.c
48 48

  
49 49
bin_PROGRAMS = \
50 50
	do_dssp		editconf	eneconv		\
......
72 72
	g_tcaf      	g_traj      	g_tune_pme   \
73 73
	g_vanhove	g_velacc    	\
74 74
	g_clustsize 	g_mdmat     	g_wham		g_kinetics	\
75
	g_sigeps
75
	g_sigeps	g_correl
76 76

  
77 77

  
78 78
LDADD = $(lib_LTLIBRARIES) ../mdlib/libmd@LIBSUFFIX@.la \
src/tools/g_correl.c
1
/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
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 4.0.99
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-2008, 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
 * Groningen Machine for Chemical Simulation
34
 */
35

  
36
/*
37
 * g_correl.c (c) 2010 Alexey Shvetsov
38
 */
39

  
40
#ifdef HAVE_CONFIG_H
41
#include <config.h>
42
#endif
43

  
44
#include <gmx_ana.h>
45

  
46

  
47

  
48
/*
49
 * This is just a wrapper binary.
50
 */
51

  
52
int main(int argc, char *argv[]) {
53
  gmx_correl(argc,argv);
54
  return 0;
55
}
56

  
src/tools/gmx_correl.c
1
/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
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 4.0.99
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-2008, 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
 * Groningen Machine for Chemical Simulation
34
 */
35

  
36
/*
37
 * gmx_correl.c (c) 2010 Alexey Shvetsov 
38
 */
39

  
40
#ifdef HAVE_CONFIG_H
41
#include <config.h>
42
#endif
43
#include <math.h>
44
#include <string.h>
45
#include <time.h>
46

  
47
#if ((defined WIN32 || defined _WIN32 || defined WIN64 || defined _WIN64) && !defined __CYGWIN__ && !defined __CYGWIN32__)
48
#include <direct.h>
49
#include <io.h>
50
#endif
51

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

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

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

  
138
  CopyRight(stderr,argv[0]); 
139
  parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
140
		    NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv); 
141

  
142
  clear_mat(zerobox);
143

  
144
  fitfile    = ftp2fn(efTPS,NFILE,fnm);
145
  trxfile    = ftp2fn(efTRX,NFILE,fnm);
146
  ndxfile    = ftp2fn_null(efNDX,NFILE,fnm);
147
  averfilea  = ftp2fn(efSTO,NFILE,fnm);
148
  averfileb  = ftp2fn(efSTO,NFILE,fnm);
149
  asciifile  = opt2fn_null("-ascii",NFILE,fnm);
150
  xpmfile    = opt2fn_null("-xpm",NFILE,fnm);
151

  
152
  read_tps_conf(fitfile,str,&top,&ePBC,&xref,NULL,box,TRUE);
153
  atoms=&top.atoms;
154

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

  
166
  if (bFit) {
167
    snew(w_rls,atoms->nr);
168
    for(i=0; (i<nfit); i++)
169
	w_rls[ifit[i]]=1.0;
170
    }
171

  
172
  /* Prepare reference frame */
173
  if (bPBC)
174
   rm_pbc(&(top.idef),ePBC,atoms->nr,box,xref,xref);
175
  if (bFit)
176
   reset_x(nfit,ifit,atoms->nr,NULL,xref,w_rls);
177

  
178
  snew(xa,natomsa);
179
  snew(xb,natomsb);
180
  snew(xava,natomsa);
181
  snew(xavb,natomsb);
182

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

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

  
320
  thanx(stderr);
321

  
322
  return 0;
323
}
0
-