Project

General

Profile

0001-Add-g_correl-tool.patch

Rebased patch - Alexey Shvetsov, 08/11/2010 03:32 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     
......
44 44
    g_angle g_bond g_bundle g_chi g_cluster g_confrms g_covar
45 45
    g_current g_density g_densmap g_dih g_dielectric
46 46
    g_helixorient g_principal g_dipoles g_disre g_dist
47
    g_dyndom g_enemat g_energy g_lie g_filter g_gyrate
47
    g_dyndom g_enemat g_energy g_lie g_filter g_gyrate g_correl
48 48
    g_h2order g_hbond g_helix g_mindist g_msd g_morph g_nmeig
49 49
    g_nmens g_order g_polystat g_potential g_rama g_rdf g_rms
50 50
    g_rmsf g_rotacf g_saltbr g_sas g_select g_sgangle g_sham g_sorient
src/tools/CMakeLists.txt.orig
1

  
2
add_library(gmxana 
3
            autocorr.c      expfit.c        polynomials.c   levenmar.c      
4
            anadih.c        pp2shift.c      dlist.c         
5
            eigio.c         cmat.c          
6
            eigensolver.c   nsc.c           
7
            hxprops.c       fitahx.c        
8
            geminate.c
9
            gmx_analyze.c   gmx_anaeig.c    gmx_angle.c     gmx_bond.c      
10
            gmx_bundle.c    gmx_chi.c       gmx_cluster.c   gmx_confrms.c   
11
            gmx_covar.c     gmx_current.c   
12
            gmx_density.c   gmx_densmap.c   gmx_dih.c       
13
            gmx_dielectric.c        
14
            gmx_kinetics.c  gmx_spatial.c   gmx_tune_pme.c
15
            gmx_dipoles.c   gmx_disre.c     gmx_dist.c      gmx_dyndom.c    
16
            gmx_enemat.c    gmx_energy.c    gmx_lie.c       gmx_filter.c    
17
            gmx_gyrate.c    gmx_h2order.c   gmx_hbond.c     gmx_helix.c     
18
            gmx_mindist.c   gmx_msd.c       gmx_morph.c     gmx_nmeig.c     
19
            gmx_nmens.c     gmx_order.c     gmx_principal.c 
20
            gmx_polystat.c  gmx_potential.c gmx_rama.c      
21
            gmx_rdf.c       gmx_rms.c       gmx_rmsdist.c   gmx_rmsf.c      
22
            gmx_rotacf.c    gmx_saltbr.c    gmx_sas.c       gmx_sdf.c       
23
            gmx_select.c
24
            gmx_sgangle.c   gmx_sorient.c   gmx_spol.c      gmx_tcaf.c      
25
            gmx_traj.c      gmx_velacc.c    gmx_helixorient.c 
26
            gmx_clustsize.c gmx_mdmat.c     gmx_wham.c      
27
            correl.c        gmx_sham.c      gmx_nmtraj.c    
28
            gmx_trjconv.c   gmx_trjcat.c    gmx_trjorder.c  gmx_xpm2ps.c    
29
            gmx_editconf.c  gmx_genbox.c    gmx_genion.c    gmx_genconf.c   
30
            gmx_genpr.c     gmx_eneconv.c   gmx_vanhove.c   gmx_wheel.c     
31
            addconf.c       calcpot.c       edittop.c       gmx_bar.c
32
            gmx_membed.c    )
33

  
34
target_link_libraries(gmxana md ${GMX_EXTRA_LIBRARIES} ${GSL_LIBRARIES})
35
set_target_properties(gmxana PROPERTIES OUTPUT_NAME "gmxana${GMX_BINARY_SUFFIX}")
36

  
37
# List of programs with single corresponding .c source file,
38
# used to create build rules automatically.
39
#
40
set(GMX_TOOLS_PROGRAMS
41
    do_dssp editconf eneconv genbox genconf genrestr g_nmtraj 
42
    make_ndx mk_angndx trjcat trjconv trjorder g_wheel 
43
    xpm2ps genion g_anadock make_edi g_analyze g_anaeig
44
    g_angle g_bond g_bundle g_chi g_cluster g_confrms g_covar
45
    g_current g_density g_densmap g_dih g_dielectric
46
    g_helixorient g_principal g_dipoles g_disre g_dist
47
    g_dyndom g_enemat g_energy g_lie g_filter g_gyrate
48
    g_h2order g_hbond g_helix g_mindist g_msd g_morph g_nmeig
49
    g_nmens g_order g_polystat g_potential g_rama g_rdf g_rms
50
    g_rmsf g_rotacf g_saltbr g_sas g_select g_sgangle g_sham g_sorient
51
    g_spol g_sdf g_spatial g_tcaf g_traj g_tune_pme g_vanhove
52
    g_velacc g_clustsize g_mdmat g_wham g_kinetics g_sigeps g_bar
53
    g_membed
54
    )
55

  
56

  
57

  
58
foreach(TOOL ${GMX_TOOLS_PROGRAMS})
59
    add_executable(${TOOL} ${TOOL}.c)
60
    target_link_libraries(${TOOL} gmxana)
61
    set_target_properties(${TOOL} PROPERTIES OUTPUT_NAME "${TOOL}${GMX_BINARY_SUFFIX}")
62
endforeach(TOOL ${GMX_TOOLS_PROGRAMS}) 
63

  
64

  
65
install(TARGETS ${GMX_TOOLS_PROGRAMS}
66
	gmxana DESTINATION ${LIB_INSTALL_DIR}	
67
      	RUNTIME DESTINATION ${BIN_INSTALL_DIR})
68

  
69
configure_file(${CMAKE_CURRENT_SOURCE_DIR}/libgmxana.pc.cmakein ${CMAKE_CURRENT_BINARY_DIR}/libgmxana.pc @ONLY)
70
install(FILES ${CMAKE_CURRENT_BINARY_DIR}/libgmxana.pc DESTINATION ${LIB_INSTALL_DIR}/pkgconfig)
src/tools/Makefile.am
46 46
	gmx_editconf.c	gmx_genbox.c	gmx_genion.c	gmx_genconf.c	\
47 47
	gmx_genpr.c	gmx_eneconv.c	gmx_vanhove.c	gmx_wheel.c	\
48 48
	addconf.c 	addconf.h	gmx_tune_pme.c  gmx_membed.c    \
49
	calcpot.c 	calcpot.h 	edittop.c
49
	calcpot.c 	calcpot.h 	edittop.c	gmx_correl.c
50 50

  
51 51
bin_PROGRAMS = \
52 52
	do_dssp		editconf	eneconv		\
......
74 74
	g_tcaf      	g_traj      	g_tune_pme   \
75 75
	g_vanhove	g_velacc    	g_membed      \
76 76
	g_clustsize 	g_mdmat     	g_wham		g_kinetics	\
77
	g_sigeps
77
	g_sigeps	g_correl
78 78

  
79 79

  
80 80
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
  gmx_rmpbc_t  gpbc=NULL;
127

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

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

  
143
  clear_mat(zerobox);
144

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

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

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

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

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

  
180
  snew(xa,natomsa);
181
  snew(xb,natomsb);
182
  snew(xava,natomsa);
183
  snew(xavb,natomsb);
184

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

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

  
323
  thanx(stderr);
324

  
325
  return 0;
326
}
0
-