|
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 |
|
-
|