0001-Add-g_correl-tool.patch
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 |
- |